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Chapter  1 


INTRODUCTION 

The  steady  operation  of  a  liquid-propellant  rocket 
engine  Is  often  disturbed  by  the  occurrence  of  large 
pressure  oscillations  in  the  combustion  chamber.  These 
oscillations,  which  can  lead  to  damage  to  or  failure  of  the 
motor,  are  caused  by  ampllficat icn  of  initially  small 
acoustic  disturbances  due  to  the  energy  released  by  unsteady 
burning  of  the  propellant.  This  is  usually  referred  to  as 
combustion  instability.  It  Is  generally  accepted  that 
unsteady  burning  can  be  correlated  with  the  gas  pressure 
(pressure  sensitivity)  and  the  magnitude  of  velocity  of  the 
fuel  drops  relative  to  the  gas  (velocity  sensitivity). 

This  dissertation  is  concerned  with  mathematical  modeling 
of  velocity-sensitive  combustion  instability  in  liquid- 
propellant  rocket  motors. 

A  survey  of  literature  dealing  with  mathematical 
modeling  of  pressure-sensitive  combustion  Instability  can 
be  found  in  the  dissertation  of  Powell  (1).  One  of  the 
first  papers  to  examine  nonlinear  effects  is  that  of  Maslen 
and  Moore  (2)  who  considered  only  the  fluid  mechanical 
effects  in  a  circular  cylinder.  The  paper  by  Priem  and 
Guentert  (3)  is  one  of  the  first  to  include  the  effect  of 
velocity  sensitivity.  They  discussed  combustion  instability 
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2 


in  a  thin  annulus. 

Since  direct  numerical  solution  of  the  governing 
equations  is  very  time  consuming  and  difficult  to  extract 
information  from,  many  investigators  such  as  Powell  and 
Ginn  (1),  Lores  and  Zinn  (5),  and  Pedciesion,  Ventrice,  and 
Purdy  (6)  have  employed  methods  of  weighted  residuals  such 
-is  Galerkin  ani  collocation  methods. 

Previous  applications  of  the  method  of  weighted 
residuals  to  combustion-instability  problems  have  dealt  with 
pressure-sensitive  combustion.  In  this  work,  this  method  will 
be  extended  to  handle  situations  in  which  velocity  sensitivity 
is  important.  Both  the  collocation  and  Galerkin  methods  will 
be  considered.  Stability  boundaries  and  pressure  wave  forms 
will  be  computed  numerically  for  transverse  motion  in  a 
cylindrical  combustion  chamber.  In  addition  a  finite- 
difference  method  will  be  used  to  solve  the  one-dimensional 
problem,  of  transverse  motion  in  a  thin  annular  chamber. 

These  results  will  be  compared  with  those  obtained  by  a 
method  of  weighted  residuals  solution  of  the  same  problem 
in  order  to  assess  the  accuracy  of  the  approximate  solution. 

In  this  study,  attempts  to  solve  the  velocity- 
sensitive  combustion  instability  problem  by  various 
numerical  methods  are  considered.  The  governing  equations 
that  describe  the  flow  conditions  inside  liquid-propellant 
rocket  motors  v/ill  be  derived  in  Chapter  Two.  An  analytical 
solution  is  obtained  in  Chapter  Three  for  nonlinear  acoustic 
motion  in  a  cylindrical  chamber.  The  collocation  method  is 
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Chapter  2 


COMBUSTOR  EQUATIONS 

*.  is  the  objective  of  the  discussion  presented  in 
pter  to  analyze  the  combustor  conservation  equation 
,:e  them  to  a  tractable  system.  The  simplification 
i  ^r.e  in  such  a  manner  that  will  allow  the  resulting 
.s  to  retain  both  the  mathematical  and  physical 
-f  the  original  problem. 

:r.der  the  assumptions  that  the  usual  balance 
_es  of  mechanics  car.  Le  applies  separately  - r  eacr. 
t he  followin'*  equations  are  derived  by  applying  the 
or.s-.-rvat ion  of  mass,  momentum,  and  energy  to  an 
ry  stationary  control  volume. 

Conservation  of  mass  of  the  fluid  requires  that  the 

* 

to  of  change  of  mass  in  a  volume  v  equals  the  rate 
r.  mass  enters  v  plus  the  rate  at  which  mass  is 
ed  inside  v  and  can  be  expressed  as  follows  (for  a 
olume )  : 

<>,_*/  pdv  =  -  /  pu*nds  +  /  wuv 


is  an  outward  unit  vector  normal  to  s,  p  is  the 

is  the  velocity  of  fluid,  w  is  the  rai 


n  x  o  Lj 


enerated  per  unit  volume ,  and  t  is  the 


4 


r 


’Jsir.g  the  divergence  the  cm  an  J  *  he  arsitrariness  of 
the  stationary  control  volume,  y lei  Is 


„  *  * 

*  »,*-►.  * 

3  _  #  p  +  v • ( pu )  =  w . 


(2.1) 


:imilirly,  the  balance  law  of  mass  for  fuel  phase  1; 


*  * 

*  -£  .  *  -V  .  * 

*sT  +  ?•(.:.  uT  )  =  -  w 


(2.2) 


"her-.-  s.  Is  the  density  of  fuel  and  uT  is  the  velocity  o: 


'he  bal 


iiance  law  or  linear*  momentum  states  that  the 


rate  of  change  of  linear  momentum  In  the  volume  v  equals 
t.-.e  rate  of  entry  of  linear  momentum  across  the  surface  s 
rlus  the  sum  of  all  external  forces  acting  on  the  volume  v 
and  can  be  expressed  as 


*  *  *  #  * 

*  -  .  *  v,+  *  ** 

h  */  sudv  =  -  f  pu(u-n)ds  +  /  rdv  +  /  wuTdv  -  /  Pnl: 

-  V  s  V  V  b  S 


•where  is  the  force  per  unit  volume  applied  to  the  gas  by 

* 

the  fuel  and  P  is  the  pressure  of  the  gas. 

Using  the  same  argument,  the  equation  becomes 


*  »* 


3  #(pu)  +  v-(puu)  =  P  +  wuT  -  7 P . 

^  Li 


(2.3) 


Similarly,  the  balance  of  momentum  for  the  fuel  phase  i: 


I 

I 


*  *  *  *  *  * 

3  „  (p\u)  +  ^.(pTS  O  =  -  ?  -  wu. 

t  *  L  L  L  L  U  U 


Substituting  (2.1)  in  (2.3)  yields 

*  #  *#  *  *  * 

#  _  #  * 

p(3  ,u  +  u  •  vu )  =  F  +  w(uT  -  u)  -  ^P. 
t  ^ 


( 

( 2  .  '4 ) 


( ?  .  5  ) 


Similarly  for  the  fuel  phase 

*  *  ** 

*L(3tA  +  Ul-VUl)  =  - 

The  law  of  conservation  of  energy  states  that  the 
rate  of  change  of  energy  in  a  volume  v  equals  the  rate  at 
•which  energy  enters  the  volume  v  plus  the  rate  at  which 
energy  is  generated  internally  plus  the  rate  at  which  work 
is  done  by  external  forces  and  gives 


# 

-> 

F , 


(2.6) 


*  * 

a  /  P  (  e 

*  V 


+  gu2)dv 


*  * 

;sp(e 


+  hu  2)(u*n)d? 


/ 


* 

(Pn) • uds 


+ 


*  * 

J.  _v  *  # 

f  F*urdv  +  /  w(eT 
v  h  v  h 


* 


^u£)d 


lv  + 


/ 


*  *  *  * 

■where  e  +  t;u2  is  the  energy  of  gas  per  unit  volume,  +  ^u 

* 

is  the  energy  of  fuel  per  unit  volume,  and  Q  is  the  heat 
transfer  rate  from  the  fuel  to  the  gas. 

Following  the  same  argument,  the  equation  becomes 


fro 


7 


*  * 

+  *&u- ) )  t  v  •  ( p  *. <.•  +  SsU-Ju)  - 

t 


-V • (Pu) 


"  ;  *  *  *  .  * 

+  P.Uj  +  w(0j  +  Vu.  )  t  Q. 


(P.6) 


: lmilar ly ,  t  hr  balance  of  energy  for  the  fuel  phase  becomes 


ft  *  #  *  *  *  *  ^ 

;t  (e  +  W  ))  +  (p.  Ce  +  W  )u  ) 

t.  *  Li  L  L  ii  ^  ^ 


*  * 


=  -  £.U,  -  w ( e.  +  1,U-')  -  Q. 


(2.7) 


Pubet  Itutln,;  (2.1).  (2.2)  into  (2.6),  (2.7)  respectively 


l  el 2o 


»  # 


!  * 


p  ( .t  „(e  +  ku2)  +  u  •  V  ( e  +  ku- )) 
t 


*  « 

-v  *  ► 

V  •  (Pu) 


*  *  *  *  * 

+  ?*u  +  w (  ( e ^  +  ljUj’ ) 


#  *  * 

(e  +  hu-))  +  Q 


(2.2) 


*  # 

»  .  *  ,  #  . ,  -»  ,* 
p LC'\  *  °L  +  +  uT,W(eL 


#  *  *  * 

+  } ) )  =  -  P.UT  -  Q. 

I »  Li 


( 2 . 0  ) 


11  la  more  convenient  to  work  with  the  thermal  energy 
equation  for  each  phase.  For  gas  phase  this  Is  done  by 


dotting  u  into  (2.3)  to  obtain  the  mechanical  energy 


equation  and  subtracting  this  from  (2.8)  and  using  (2.1)  to 
get 


*r.  *  # 

pDe/Dt  — 


PD  /Dt 


*(logp) 


« 


(  Ur  -  u)  +  Q  + 


*, ,  * 

w((eT  - 


*  > 

e) 


*  2 

(u.  - 


*  2 

u  )/2 


V  ^ 

+  u • (u  -  UT  )  - 


*  * 

P/d) 


(2.10) 


In  a  similar  way,  it  can  be  shown  that 


* 

P 


D  e  /Dt 
L  L  L 


* 


* 

Q 


(2.11) 


where  D  /Dt* 


*  * 

u  •  y , 


dt  /Dt* 


*  # 

3  +  u  •  ^  are  the 

t*  L 


comoving  derivatives  for  each  phase. 

For  simplicity,  gas  phase  viscosity  and  heat 
conduction  were  neglected  in  the  previous  analysis.  Since 
these  produce  dissipation,  equations  derived  in  this  way 
should  give  a  conservative  estimate  of  the  stability  of  the 
system. 

It  appears  (see,  for  instance,  Powell  (1)}  that  the 
primary  effect  of  the  burning  fuel  is  that  of  interphase 
mass  transfer.  Thus, the  balance  laws  for  mass,  momentum, 
and  thermal  energy  for  each  phase  will  be  further 
simplified  by  neglecting  all  interphase  transfer  terms 
other  than  those  appearing  in  the  mass-balance  equations. 


This  leads  to  the  system  of  equations: 


Dp/Dt*  + 


,  *  * 


* 

w 


»  # 

D*r  /Dt  +  Pj  V  -u.  =  -  w 


*  ** 

*  ►  ,  *  *r, 
pDu/Dt  =  -  VP 


*  ►  ,  * 

PrDTuT/Dt  - 


'L  L  L 


0 


pCvDT/Dt*  =  P(D  /Pt*)(logp) 

*  *  * 

p  C  DT  /Dt  =  0  (2.12) 

L  vb  r, 

*  *  *  * 

where  the  constitutive  equations  e  =  C  T  and  e^  =  C^Tj 

have  been  used  (with  T  denoting  temperature  and  (  ^  spec  ■  *  1c 
heat).  The  equation  of  state  for  an  ideal  gas  is 


* 

P  = 


* 

pRT 


(2.13) 


where  R  Is  the  gas  constant. 

The  governing  equations  will  now  be 
with  respect  to  a  steady  reference  state, 
denoted  by  the  subscript  "r" .  Ail  lengths 
to  some  characteristic  length  i-r,  such  as 


no  nd i mens Penalized 
which  will  be 
will  be  referred 
the  chamber  radius. 


The  characteristic  velocity  is  the  speed  o!’  sound  at  the 
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re; erence  state,  and  the  characteristic  time  is  the  wave 
travel  time  Lr/Cr>  The  dimensionless  quantities  are 
defined  below: 


V  =7  /L 


i  > 


L^t/C 


u  =  Cru 


UL=  CA 


o  =  orp 


PL=  °r°L 


p  =  p  p 

rr 


w  =  orCrw/Lr 


T  =  C„T/(yR) 


TL=  CrTL/(YR) 


Y?r>/  P> 


^r»  Pv*C>,/y 


rhe  dimensionless  conservation  equations  become 


Cont inui 


uIZ 


Dp/Dt  +  p  V • u  = 


Dpr/Dt  +  P  ^.u  =  -  W 

i-j  L 


(2.14) 


Momentum 


Du/Dt  =  -  ^P/Y 
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DuT/Dt  =  0 


(2.15) 


.".nergy 


DT/Dt  =  (y  -  1 )  P  (  D  /Pt  )  (logp) 


DTj/Dt  =  0 


(2.16) 


Equation  of  state 


F  =  PT 


(2.17) 


By  solving  (2.16a)  one  obtains 


T  =  p>-1 


(2.18) 


Equation  (2.17)  then  become; 


P  *  p 


(2.19) 


Substituting  (2.19)  to  (2.15a)  produces 


Du/Dt  =  -v(p^-1)/(y-  1)  •  (2.20) 

Taking  the  curl  on  both  sides  of  equation  (2.20)  gives 

$  X  Du/Dt  =  0  .  (2.21) 

It  can  be  shown  using  vector  identity  that  (2.21) 
leads  to  the  following  equation  which  describes  the 
generation  of  the  vortlcity  0  =  $  X  u  in  the  flow. 


Dfi/Dt  =  Q  •  (  V  u )  -  ( v  •  u  )  . 


(2.22) 
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This  equation  describes  the  variation  of  vorticity 
for  a  given  fluid  particle.  Suppose  that  at  t  =  0 ,  this 
particle  has  zero  vorticity  and  at  the  point  where  it  is 
located  the  velocity  gradients  are  bounded.  Then  by  (2.22) 
the  rate  of  change  of  vorticity  of  the  particle  vanishes. 

It  follows  from  this  that  the  vorticity  will  be  zero  at  the 
next  instant.  3y  induction,  it  is  seen  that  as  long  as  the 
velocity  gradients  are  bounded  at  each  point  occupied  by  the 
particle,  its  vorticity  will  remain  zero  for  all  time.  If 
ail  fluid  particles  in  the  system  have  zero  vorticity  at 
t  =  0,  the  vorticity  will  vanish  at  all  points  in  the  flow 
field  and  for  all  t  >_0,  therefore,  as  long  as  the  velocity 
gradients  are  bounded  through  the  flow  field  (  as  they  would 
not  be,  for  example,  at  a  shock  wave  ).  Assuming  this 
condition  to  be  satisfied  implies  irrotational  flow  (n  =  0) 
and  allows  the  definition  of  a  velocity  potential  <j>  defined 
such  that 

U  =  V*  .  (2.23) 

Substituting  (2.23)  to  (2.20)  gives 

3  $  +  4>  -  v  4>  +pl~1/(y-l)  =  F(t)  . 

Since  the  velocity  is  the  space  derivative  of  <j> ,  it 
is  permissible  to  add  to  <P  any  arbitrary  function  of  time. 

This  is  equivalent  to  a  statement  that  F(t)  is  arbitrary. 

For  future  convenience  F(t)  was  selected  to  be  l/(y-l). 

This  choice  yields 

O 


oT  J=  1  -  C y— i ) C at ♦  +  Hv<> •  v 4> j 


Substituting  (2.2^)  into  (2.14a),  (2.18)  and  (2. in)  gives 


3  $  -  V  2  ^>  (  1  -  (  y  - 1 )  (3  $  +  H’v<>  •  V<p )  )  +  V  $  •  V  (  2  3 

t  t  t,  t 


+  'i  V  $  •  V $ )  +  (1  —  (  y  — l )  (  3  4>  +  ^Vi^,Vcj)’'  ^  ^  ^  ^  W  -  0 


T  =  1  -  (y— l)(a,4>  +  $• $$) 


?  =  (1  -  (  Y-l )  (  3.  ♦  +  SV^-^))Y  ,(y_:  ' 


This  system  of  equations  contains  nonlinear  .-as 
dynamics  terms  of  all  orders. 

Due  to  their  highly  nonlinear  and  mathematically 
complicated  nature,  the  system  of  equations  obtained  abov 
cannot  be  solved  exactly.  In  order  to  obtain  simpler,  bu 
approximate,  equations  which  can  be  more  easily  analysed, 
nonlinear  terms  of  order  higher  than  two  are  neglected. 
This  has  the  effect  of  retaining  the  most  important 
nonlinear  effect,  while  greatly  simplifying  the  form  of  t 
equations.  The  system  of  equations  then  becomes 


1  -  (  3 1 >(>  +  ‘,v*.y\.-.>)  +  \,(2-y)  (  3t$)2 


T  =  1  -  (y-1)  (  3.  4>  + 


P  =  1  -y(at?  +  <J>  •  V  <|>  )  +  \y  (  9  $  )  2 


1 4>-V24,  (1  -  (y-l)3  <>)  +  270.7(3  <>) 


+  (1  -  ( y-2)  3  >)«  =  0 


(2.26) 


For  a  cylindrical  combustor,  it  is  desired  to 
investigate  the  stability  of  the  steady-state  solution  of 
(2.26d).  To  do  this  the  steady-state  solution  must  be 
obtained.  It  is  assumed  that  the  steady-state  burning  rate 
depends  on  z  only  (w  =  w(z))  and  that  the  flow  is  axial 
($  =  £_(z)).  Thus  (2.26)  simplifies  to 


du/d  z  =  w 


(u  =  d£/d  z) 


T  =  1  -  ?s(Y-l)u2 


P  =  1  -  byu2 


=  1  -  Jsu2  . 


(2.27) 


It  can  be  seen  from  (2.27)  that  the  deviations  of  the 
steady-state  solution  from  a  uniform  state  are  0(u2).  To 
allow  a  discussion  of  orders  of  magnitude  of  various  terms 
define  the  ordering  parmeter  e  such  that 


e  =  max  (u) . 


(2.28) 


lk 
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his  Imp lie 


;  hat 


u  =  0(e),  w  =  0(c),  J  =  0(c).  (2. 29) 

The  stability  analysis  will  now  be  carried  out  by  assuming 
that 


=  ±(  z)  +  '  (  x,  y,  z,t ) 

W  =  w( 2 )  +  e  2w  ' ( x, y, s,t ) 

(2.30) 

w’  =  0(1). 

(2.21) 

It  is  being  assumed  that  the  unsteady  perturbation  of  the 
velocity  potential  from  the  steady  state  is  of  the  same 
order  of  magnitude  as  the  deviation  of  the  steady  state 
from  a  uniform  state  and  that  the  unsteady  burning  rate 
perturbation  is  of  the  sane  order  as  the  unsteady  nonlinear 
gas-dynamic  terms.  The  first  assumption  is  necessary  to 
be  consistent  with  the  quadratic  approximation  inherent  in 
(2.26)  while  the  second  eliminates  nonlinear  terms  involving 
products  of  w'  and  These  assumptions  are  made  for 

simplicity  and  are  not  meant  to  imply  that  other  orders  of 
magnitude  for  the  various  terms  are  impossible  or  even 
unlikely.  The  objective  of  this  work  is  to  pick  one  case 
which  is  mathematically  tractable  and  subject  it  to  extensive 
analysis.  While  it  is  felt  that  most  of  the  features  of  the 
stability  problem  will  be  reflected  by  this  case,  it  is 
recognised  that  many  other  possibilities  remain  to  be 
explored . 
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Substituting  (2.30)  and 

p  =  £  +  eo',  P  =  P  +  Ep ' ,  T  =  T  +  eT'  (2.32) 

(o'  =  0(1),  P»  =  0(1),  T'  =  0(1)}  into  (2.26), retaining  all 
terms  of  Q(e3)  or  lower, and  dividing  all  resulting  equations 
by  e  leads  to 

3tt»'  -  v2*'  +  ^st*'  + 

+  c{(y- 1 )  V  2  0  »  3t  4>  ’  +  2^’  *V(  3t*' )  +  w»>  =  0 

o'  =  -  {3t^'  +  U3^'  +  te(vV  *V$'  -  (2-y  )  (  3>t  if  '  )  2  )  > 

P'  =  -  t(3j'  +  u3  $’  +  \>  c  (  V  $  1  *  V  4> '  -  (3  *')2)} 
t  ~  z  t 

T'  =  -  (  y-1  )  (  3 1  <»> '  +  U3J'  +  VeV^'-V^').  (2.33) 

The  equations  derived  from  this  perturbation  analysis 
are  expected  to  be  valid  as  long  as  the  amplitudes  of  the 
perturbation  quantities  are  finite  and  smaller  than  unity. 

For  simplicity  the  primes  appearing  in  the  above 
equations  will  be  dropped  and  it  will  be  understood  that 
unprimed  quantities  are  associated  with  perturbation  from 


he  steady  state. 


Chapter  3 


ANALYTICAL  SOLUTION  OF  WAVE  EQUATION 

In  this  chapter  a  second-order  solution  for  the 
velocity  potential  associated  with  transverse  acoustic  wave 
notion  of  an  ideal  gas  in  a  circular  cylinder  is  obtained. 
The  solution  is  found  in  the  form  of  a  Fourier-Bessel 
series.  This  provides  a  standard  to  compare  the  proposed 
methods  which  will  be  discussed  in  the  following  chapters. 

The  governing  equation  is  obtained  from  (2.33a)  by 
assumption  that  the  burning  rate  function  Is  absent.  Then 
the  equation  reduces  to  the  second-order  form  of  the  one 
discussed  by  Maslen  and  Moore  (2)  for  pure  gas  motion  in  a 
cylinder. 


-  V 2  <i>  +  e(2v<j)  •  V  3t  <t>  +  (y-1)v24i  3t4>)  =  0  (3.1) 

For  a  cylinder  of  radius  Lr,  a  set  of  cylindrical 
polar  coordinates  (r,o,z)  with  the  z  axis  being  coincident 
with  the  cylinder’s  axis  of  symmetry  is  used.  In  seeking 
a  second-order  transverse  wave  solution,  one  assumes 

4>  =  *  ( r ,  0  ,  t )  +  e  j>  ( r  ,  e  ,  t )  +  ...  .  (  3  •  2 ) 

l  2 

Substituting  (3-2)  into  (3.1)  and  (2.33)  loads  to 


' » 

/ 
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tt. 


$  -  V2?  =  0 


a-t*2  -  '}2'i 


-  (a,(u2)  +  (y-l)v2*  M  )  (3.3) 

c  i  l  z  l 


P  =  1  +  CP  +  E  2  P 
1  2 


(3.4) 


where  ?,  =  -  y3  ?  ,  P  =  +  ^(u2  -  (  3  <?  )2))  (3*5) 

1  t  l  2  z  2  l  r  l 


72>i  =  3rr*i  +  3r*i/r  +  396*i/r2’  1  =  X*  2 


u2  =  (  ^  ) 2  +  (jg^/r)2. 


(3.6) 


Three  different  cases  are  discussed  below.  For  the 
initial  conditions 


}(r,9,0)  =  J  (S  r ) cos 9  ,  3.$(r,e,0)  =  0 

ill  1 


(3.7) 


the  solution  of  (3- 3a)  is 


a  =  J  (3  r)cosecos(S  t), 
l  in  ll 


Substituting  this  equation  to  (3- 3b)  produces 


3..*  -72?  = (b  +  i  l  Jm(Smnr)cos (m9 ) )sin(2S  t)  (3-8) 
tt  *  2  2  oo  m=  o  n=  i  111  ll 

l 

•where  J  is  used  to  denote  an  n'th  order  Bessel  function  of 
n 

the  first  kind,  Smn  is  used  to  denote  the  m’th  root  of 
the  equation  =  0 ,  and  a  prime  indicates  differentiation 


N 


.  — .. 


**  T 


With  re3?ect  to  the  argument  of  the  Ber.sel  :  met  Ion.  The 

coefficients  b  are  intuTPii-  o  1  „ 

mn  aij  o.  Bessel  functions  which 

:-.ust  be  computed  numerically.  The  details  of  this 

calcuat ion  are  discussed  in  appendix  A. 

Solving  (2-3)  one  obtains 
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L 


2  =  (b00  -  sin(2Sut)) 


+m=  o  n=  i  (bmn/(3mn  '  ,  ))Jra(STinr)cos  ^  -  )  (sin(2S  t) 

1 1  “  i  i 


m*  i 


-  (23  /s  )sin(S  t)) 

1 1  mn  mn  '  • 


(3.9) 


"e  °““C’r  tv/°  £C*utions  to  be  discussed  subsequently  are 
-"ur.J  In  the  same  way.  For  the  initial  conditions 

*(r’  s»°)  =  °>  \4(r,Q,0)  =  s  J  (S  r)cosa  (^.13) 

L  iiiu  u 

t  is  found  that 


=  J : 1r)cos0sin(3  t) 


(3-11) 


ind  tnat  ^ 2  is  eiven  by  the  negative  of  (3.9).  For  the 
initial  conditions 


*(r,9,0)  =  J  (3  r)cose, 

ill  ’ 


^♦(r,e,0)  =  3  J  (S  r)sinO 

ill  li 


(3-12) 


;he  solution  Is 


_ ' 


J 
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$  =  J^CS^^rJcosCo-i-it  -  G  ) 


IV 


and  ^ 2  =  rJi{2b2n/(S2n  }J2  (S2nr  >  (sin( 26  >  (cos  (2Snt ) 


-  cos(S2nt)}  +  cos(29)  {sin(2S11t)  -  (2S11/S2n)Sin(S2nt ) } ) 


It  can  be  seen  that  for  these  initial  conditions  $ 
represents  a  spinning  wave  but  4  j  +  e<p2  does  not.  Some 
typical  results  for  the  pressure  functions  P,  and  P_  were 

i-  d 

computed  using  the  solutions  for  initial  conditions  (3.7) 
(labeled  standing  wave)  and  initial  conditions  (3.12) 

(labeled  traveling  wave)  for  r  =  1,  6  =  o,  and  y  =  1.2. 

These  data  are  presented  graphically  in  Figure  1.  For 
these  computations  the  series  expansions  were  terminated  at 
n  =  5.  There  is  virtually  no  difference  between  the  results 
for  n  =  4  and  n  =  5,  and  even  the  results  for  n  =  1  provide 
a  reasonably  accurate  solution. 

The  results  discussed  above  were  used  to  check  the 
numerical  calculations  to  be  discussed  in  the  next  two  chap¬ 
ters.  These  results  generalize  those  of  Maslen  and  Moore  (2) 
to  initial  conditions  which  do  not  lead  to  periodic  solutions. 
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Chapter  4 


-< 


COLLOCATION 

In  previous  investigations  of  velocity-sensitive 
combustion  instability  (see,  for  instance,  (6)  and  the 
references  therein)  the  most  widely  used  burning-rate 
equation  is  that  associated  with  the  vaporization  limited 
combustion  model  discussed  by  Priern  and  Ouentert  (3)*  In 
the  notation  of  the  present  thesis  the  second-order  version 
of  this  equation  can  be  expressed  as 


w 

= 

(w/ 

£•)((!  -  ht. 

V 

)((1 

-  u  /u 

ry 

L>’+(ut 

'uL)2)  ' 

-1) 

(4.1) 
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1 1  u 

do  of  the 

droplet 

ve 

loc 

ity 

vector  which 
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to  be 

a  x  i  a 

I  t 

Sub st  i 

tut 

ing 

(4 

.  1) 

into 

(2.33a)  yie 

Ids 

3tt*  +  St  ®t,  ♦  + 

2u 

$ 

-  v?*(l  - 

t  (  Y 

-1)  ) 

+  3  4  +  (w/f  }((1  -  \.e  'A  if )  ( ( 1  -  5  ^/u,  )  “ 

b  v  2 

+  ( v  *  v  «>  -  (  )  •'  )/u.2  1)  =  0.  ( 4 .  o  ) 

I’  1  j 
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No  exact  solution:-,  are  to  be  expected  to  (•'<.?).  Fully 
numerical  solutions,  on  the  other  hand,  are  expectec  to  be 
quite  time  consuming  for  multidimensional  problems.  For 
these  reasons  it  war.  decided  to  attempt  approximate 
partially  analytical  solutions  using  the  method  of  weighted 
residuals  (see  for  instance,  (6)1.  For  problems  of  pressure 
sensitive  combustion  instability  Powell  (1)  successfully 
employe!  the  Galerkln  method.  This  method  !s  limited  in 
•.pplicab  llit.y ,  however,  to  equations  containing  nonlinear- 
itier  involving  only  polynomial  functions  of  the  dependent 
variable  and  its.  derivatives.  It  is  clear  that  (JJ.P)  is  not 
of  this  form.  The  orthogonal  collocation  method  was, 
therefore,  chosen  because  of  its  simplicity  and  adaptability 
to  equations  containing  nonlinearities  of  any  algebraic 
form..  The  application  of  this  method  to  the  problem  of 
transverse  wave  motion  in  a  cylindrical  chamber  will  be 
iiscussed  below. 

Consider  the  transverse  wave  motion  in  a  cylindrical 
combustor.  It  is  convenient  to  describe  this  problem  in 
terms  of  cylindrical  polar  coordinates  (r,0,s).  Then 


*  =  <J>  ( r  ,  o  ,t )  . 


('1.3) 


It  will  also  be  assumed  that  w  and  are  constants.  Then 
(3-P)  becomes 


'  p  +  wo  p  -  (  <*>  1)  +  3  $/r  +  3  j/r2)(]  -  .  (y-1)3  p) 

t.t,  -  t.  rr  r  0(1  1 


+  2r(3^$3t<>  +  So  ^  30t  ^/'1'2  ^  +  “  V  e  3^  $  )  (  (1 
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+  ((3r>)2  +  (30«/r)2)/u2)V  -  l)  =  0  .  (4.4) 


Now,  assume  an  approximate  solution  of  the  form 


P  Q 

::  z 

m= on=  i 


Jm(; 


;  r) cos (mo ) f  (t) 
mn  ran 


( 4 .  5 ) 


Equation  (4.5)  is  a  superposition  of  the  normal  modes 
associated  with  the  corresponding  linear  acoustic  problem. 

A  solution  of  this  form  allows  one  to  investigate  the 
Influence  of  nonlinearities  on  the  behavior  of  modal 
amplitudes  which  would  exhibit  simple  harmonic  motion  in 
the  linear  case.  Only  standing  modes  can  be  described  by 
(4.5).  To  represent  spinning  modes  another  series  involving 
sln(me)  must  be  added.  This  is  omitted  for  simplicity  in 
this  discussion.  It  is  convenient  to  rewrite  (4.5)  as 


N 

T.  J 

J  =  l  "‘j 


r ) cos (m . 6 ) f . (t ) 
<  I  J 


(4.6) 


where  N  =  (P  +  1)Q  . 

The  equation  is  now  evaluated  at  N  points  (the  collocation 
points)  to  yield 


4> 


i 


CiJfJ 


1 


,2,.  .  .N 


(4.7) 


where 


(4.8) 


Cij  =  Jm  (Sm  njri)cos(mJ®i)' 


Next,  (4.7)  Is  inverted  to  obtain 


N  , 

h  V ,hcJv 


(4.9) 


Then  the  various  derivatives  of  (4.6)  can  be  expressed  in 
terns  of  value  of  $  as 


(V>i  “  j.VuV 


’  <3oo',)i  ~  .Jhu’j  o.io 


i  here 


r 

=  (3  C„ 

)  c  1 

ij 

Vk 

kj 

0 

=  (30Cik 

i 

,  -1 

ij 

)Ckj 

‘  i1  1 

3  C  ) C~ 1 
0,0?  ik  kj 


(4.11) 


Substituting  into  (4.4),  a  set  of  N  ordinary  second-order 
differential  equations  having  the  following  form  is 
generated . 


i 


N  .  N  N 

4,  +  w$.  -  I  C  .(I-c(y-I)*.  )  +  2e  J  E  C,  .  d>  a 
1  j  =  i  i-J  J  J  j  =  i  k=i  xJk  J  k 

N  N  t. 

+  (w/c )  (  (l-he<K  )  (1  +  Z  Y.  C  *  ♦  /u2)  --1)«0  (4.12) 

1  j  =  i  k= l  ljk  J  k  L 


where 


ij  ~  CIJ  Cij/ri  C 1  j  /ri 


■ 


25 


Cijk 


cUcik  +  c° 


PU  /r,2 

ljclk/ri 


(4.13) 


Solutions  were  obtained  by  a  fourth-order  Runge-Kutta 
me;hod  to  the  set  of  equations  (4.12).  It  was  found  that 
the  results  were  very  sensitive  to  the  number  and  locations 
of  collocation  points,  especially  to  the  radial  distribution 
of  points.  It  was  found  that  even  in  cases  when  w  was  equated 
to  zero  (no  combustion)  it  was  possible  to  find  many  choices 
of  collocation  points  which  lead  to  the  computed  modal 
amplitudes  increasing  without  bound. 

In  order  to  illustrate  the  difficulties  involved  in 
a  simple  case,  consider  a  one-dimensional  problem  of 
longitudinal  wave  motion  with  y  =  1,  u  =  0  and  w  =  0.  In 
this  case  (4.2)  reduces  to 


3  a  -  3  <)>  +  2e(3<j>3  A )  =  0 

tt  zz  z  zt 


(4.14) 


First  consider  the  boundary  condition 


3^(0, t)  =  0,  ( n  ,  t )  =  o  .  (4.15) 

A  one-term  solution  satisfying  (4.15)  Is 

$  =  f  (t)cos(-iz) .  (4.16) 

Applying  the  collocation  method  with  a  collocation  point  z0 
produces  an  equation  for  $0  =  4>(z0,t)  of  the  form 


»  vl 
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i’  +  K*  +  H,(tan2(z  /2 )  4,  4.  )  =  0 .  (4.17) 

00  000 

The  Galerkin  method  (to  be  discussed  in  detail  in  the  n^xt 
chapter)  when  applied  to  the  same  problem  yields  an  equation 
for  f  of  the  form 


f  +  >sf  +  4  ff/(3n )  =  0.  ( h . 18) 

As  another  example  consider  the  boundary  conditions 

3„4>(0,t)  =  0,  3Z  4>  (  tt  ,  t )  =  0.  (4.19) 

A  solution  satisfying  (4.19)  is 

4>  =  f(t)cos(z)  (4  20  ) 

for  which  the  collocation  method  leads  to 


4>  +<(>  -2etan2z  -f>  <}>  =  0  (4.21) 

00  000 

and  the  Galerkin  method  leads  to 

f  +  f  =  0.  (4.22) 

From  the  above  two  cases,  one  can  observe  that  the 
last  term  in  each  of  the  equations  obtained  by  collocation 
can  take  on  any  value  between  0  and  “>  depending  on  the 
location  zQ  of  the  collocation  point.  No  such  difficulty 
is  encountered  when  using  the  Galerkin  method.  F.ven  if  more 
points  were  used,  the  same  behavior  was  obtained.  Therefore 


Chapter  5 


THE  GALERKIN  METHOD 

The  Galerkin  method  Is  a  special  application  of  the 
method  of  weighted  residuals  (usually  referred  to  as  MWR). 

It  has  been  extensively  used  in  the  solution  of  various 
stability  and  aeroelasticlty  problems  {see  for  instance 
Powell  (1)  )  and  proved  itself  as  a  useful  tool  for  the 
solution  of  both  linear  and  nonlinear  problems.  Although 
it  is  an  approximate  mathematical  technique,  it  has 
nevertheless  produced  results  which  were  in  excellent 
agreement  with  available  exact  solutions.  These  approximate 
solutions  are  usually  simpler  in  form  than  the  exact 
solutions  obtained  by  numerical  integration,  and  their 
quantitative  evaluation  requires  considerably  less 
computation  time.  However,  this  method  is  known  to  be 
reliable  and  applied  conveniently  only  to  equations 
involving  nonlinearities  of  a  polynomial  type.  Therefore, 
if  this  method  is  used  in  conjunction  with  the  valorization 
limited  burning  rate  function  of  (4.1)  the  problem  becomes 
intractable.  Thus, the  following  simpler  purely  phenomenolog 
ical  burning-rate  function  was  employed  for  the  case  of 
instantaneous  combustion  response. 

w  =  wnu2  (5.1) 
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where  n  is  a  constant  which  will  subsequently  be  referred  to 
as  the  interaction  coefficient.  It  can  be  determined  either 
by  comparison  of  predictions  of  the  theory  directly  with 
experiment  or  by  comparison  with  burning  rate  laws  meant  to 
apply  to  special  types  of  combustion  processes.  As  an 
example  of  the  latter  method  consider  the  vaporization- 
limited  burning  rate  law  (4.1).  As  it  stands  (4.1)  exhibits 
both  pressure  and  velocity  sensitivity.  A  purely  velocity 
sensitive  law  can  be  obtained  by  assuming  e  <<  1  to  get 

wy  =  ( w / £ 2 ) ( ( 1  +  (u/uL)2)^  -  1)  (5.2) 

where  the  v  is  used  to  denote  the  vaporization  model.  If  one 
equates  the  slopes  of  (5.1)  and  (5.2)  at  u2  =  0  and  plots 
both  functions,  the  graph  would  look  as  follows. 


It  can  be  seen  that  (5.1)  would  overestimate  the  burning 
rate.  Thus  the  upper  bound  for  n  can  be  computed.  From 
(5.1)  and  (5.2) 

d  w  =  wn,  d  w  =  (w/(4e2u2))/(1  +  (u/u  )2)3/\ 
u  2  —  u  2  v  —  L  h 

Thus  d,w(0)=wn,  d  0w  (0)  =  w/(4e2u2). 

u‘  —  u  2  V  —  L 

Equating  these  results,  one  gets 

n  =  l/(4e2u2).  (5-3) 

max  L 

The  phenomonological  law  (5.1)  can  be  related  to  other 
special  burning-rate  laws  In  a  similar  manner. 

If  it  is  desired  to  consider  history-dependent 
combustion  processes t this  can  be  done  through  the  general 
formula 

t 

w  =  nf  G(t-s)d  (u2)dg  (5.4) 

0  5 

where  G  is  memory  function  and  £  is  a  dummy  variable.  If 
G(t)  =  H(t)  (H  being  the  unit  step  function)  all  increments 
of  change  in  u2  occurring  in  the  past  are  counted  equally 
and  (5.1)  is  recovered.  If  G(t)  =  H(t)  -  H(t-x)  all 
increments  of  change  in  u2  are  counted  equally  up  to  x  units 
of  time  in  the  past  while  those  previous  to  that  time  are 
not  counted  at  all.  Substituting  this  expression  into  (5.4) 
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and  integrating  one  obtains 

w  =  wn(u2  -  u2)  (5.5) 

—  T 

where  u  =u(t-t )  and  t  is  called  the  tine  delay.  Equation 
r 

(5.5)  will  be  used  to  account  for  the  history  of  the  burning 
process  in  a  rough  way.  More  sophisticated  treatments  of 
this  phenomenon  could  be  obtained  by  employing  a  more 
realistic  function  for  G(t).  Substituting  (5.5)  and  (5-1) 
Into  (2.33a)  one  obtains 


3 .  t  -  V2*  +  2 u 3  $  +  w  3 .  4>  +  c((y-1)V2$3  4i  +  2  V  $  •  V  3  $ 
it  —  ZZ  —  Z  Z  t 

+  wn(v$*7$-,)v$  •  v  $  ))  =  0  (5.6) 

—  T  T 

where  J  =  0  for  Instantaneous  combustion  and  j  =  1  for 
history- dependent  combustion . 

The  most  general  solution  of  (5.6)  (subject  to  hard- 
wall  boundary  conditions  for  the  unsteady  variables)  can  be 
written  in  the  form  of  the  following  Fourier-Bessel  series. 


<■  -  ‘Wt)  * 


+  "  z  (f  (t)cos(me)  +  -  (t)sin(mo)  ).T  (S  r)  (5.7) 
-=in=i  mn  mn  m  mn 


Substituting  (5.7)  into  (5.6)  and  applying  the  usual 
•1  ilerkin  orthogonalizat  ion  procedure  leads  to  an  infinite 
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set  of  coupled  ordinary  differential  equations  governing  the 
functions  fmn  and  Smn-  No  solution  can  be  obtained  unless 
the  series *(5.7)  is  truncated  so  as  to  produce  a  finite  set 
of  equations.  If  the  terms  neglected  are  actually  small, an 
accurate  solution  will  result. 

In  this  thesis  attention  will  be  focused  on  initial 
disturbances  having  the  form  of  the  first  tangential  mode. 
The  simplest  finite  series  contained  in  (5.7)  capable  of 
modeling  the  effect  of  quadradic  nonlinearities  In  (5.6)  is 


\(t)JQ(Soir)  +  f  -^t  )J1(S111’  )cos0  +  f2(t  )J^(S2ir  )cos26 


+  g1(t )J1(S11r)sin6  +  g2 ^ ^ J2^S21r^ sin2e '  (5.8) 

The  Galerkin  method  then  produces  the  following  ordinary 
differential  equations  governing  these  five  variables. 

V  s-;o  *  soiro  +  '(ciVo  +  VVi  + 

*  VVs  *  E2E2>  +  !In(ClTo  +  C12(fl  +  +  C13<f2  +  E2> 


*  C12(fl,  +  +  C13<f2,  +  S2,»»  *  0 


q  +  wfq  +  sJjfj  +  e(c11f()fl  +  c  +  r^) 


+  C7(£2E1  +  W  +  S"<CK,Vl  *  C15(flf2  +  E1E2> 


The  coefficients  C.  (  i  =  l  ,2 , .  .  .  17)  are  integrals  of 
Bessel  functions  which  must  be  computed  numerically.  The 
details  o*  thus  calculation  are  discussed  in  aoper.dix  A. 

The  stability  boundaries  in  the  (n,c)  plane 
were  letermined  first.  For  a  given  value  of  w,  a  value  of 
e  was  selected  and  solutions  were  obtained  for  various  value 
of  n.  In  each  case,  if  the  modal  amplitudes  exhibited 
growth  with  time,  the  value  of  n  was  decreased  for  the  next 


run  while  if  decay  of  modal  amplitudes  was  observed 


he 


value  of  n  was  increased  accordingly.  A  systematic 
Iteration  process  was  devised  so  that  the  computer  could 
carry  out  these  calculation:',  without  analyst  Intervention. 
In  this  way, one  point  on  the  stability  boundary  was 
established.  Then  a  new  value  of  e  was  chosen  and  the 
iteration  process  was  repeated  to  determine  another  point 
on  the  stability  boundary.  This  procedure  (which  consumed 
a  considerable  amount  of  computer  tine)  was  continued  until 
enough  points  had  been  found  to  establish  the  shape  of  the 
stability  boundary.  It  should  be  pointed  out  that  the 
amplitude  of  the  IT  mode  always  initially  decreases  due  to 
the  fact  that  purely  velocity-sensitive  combustion 
instability  is  always  linearly  stable.  Thus  it  is  necessar 
to  continue  the  calculation  for  a  cons iderable  period  of 
time  to  determine  whether  a  given  set.  of  conditions 
corresponds  to  nonlinear  stability  or  instability. 

Figures  2  and  3  show  some  typical  stability 
boundaries  for  instantaneous  combustion  using  the  initial 
condit ions 


f  ( 0  )  =  1,  fn(0) 


r .,  ( n )  =  2,(0)  =  g,(0)  =  o 


respectively.  The  region  below  and  to  the  left  of  the 
stability  boundary  is  associated  with  stable  conditions, 
while  the  region  above  and  to  the  right  is  associated  with 
unstable  conditions. 

In  the  case  of  linear  acoustics,  the  first  set  of 
initial  conditions  would  lead  to  a  standing  wave  and  the 
second  set  of  initial  conditions  corresponds  to  a  traveling 
wave.  In  both  cases, it  can  be  seen  that  the  stability 
boundaries  have  roughly  the  forms  of  rectangular  hyperbolas. 
As  expected,  increasing  the  steady-state  burning  rate  reduces 
the  value  of  n  required  to  produce  instability.  This  effect 
is  somewhat  more  pronounced  for  the  traveling  waves  than  for 
the  standing  waves.  There  does  not  appear  to  be  a  distinct 
pattern  to  the  results.  Thus  for  w  =  0.2  standing  waves  are 
more  stable  than  traveling  waves  while  for  w  =  0.1  traveling 
waves  are  more  stable  than  standing  waves. 

Figures  4  and  5  show  stability  boundaries  associated 
with  initial  conditions  (5.10)  and  (5.11)  ,  respectively,  for 
history-dependent  combustion.  It  is  apparent  that  the 
influence  of  the  time  delay  parameter  on  the  results  for 
standing  waves  is  much  greater  than  on  the  results  for 
traveling  waves.  Again,  no  clear  pattern  emerges  from  the 
results.  Comparing  Figures  2  and  4,  and  3  and  5  shows  that 
instantaneous  combustion  can  be  either  more  or  less  stable 
than  history  dependent  combustion  depending  on  the  value  of 
the  time-delay  parameter.  Comparing  Figures  4  and  5,  it 
can  be  seen  that  standing  waves  can  be  either  more  or  less 
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stable  than  standing  waves  for  history  dependent  combustion. 
Figures  4  and  5  also  illustrate  the  fact  that  increasing  the 
amount  of  time  delay  can  either  increase  or  decrease  the 
stability  of  the  system.  Thus  }  for  both  standing  and 
traveling  waves, t  =  it  corresponds  to  a  more  stable  situation 
than  either  t  =  0.5ir  or  t  =  l.5n  .  The  lack  of  patterns 
exhibited  by  these  results  emphasizes  the  need  for  numerical 
methods  of  the  type  developed  during  the  present  research. 

It  appears  that  the  only  way  to  find  out  what  will  happen  in 
a  given  situation  is  to  solve  the  equations  for  that 
particular  case.  This  matter  will  be  discussed  further 
subsequent ly . 

The  pressure  perturbation  can  be  calculated  from  the 
velocity  potential  perturbation  by  substituting  equation 
(5-8)  into  (2.33c),  expanding  the  result  in  a  Fourier-Bessel 
series  and  retaining  only  the  terms  corresponding  to  the  IT, 
2T  and  1R  modes  to  obtain 


p  = 


Y('(‘d01f0  +  L^d02f0  +  d  0  3  ^ 1  1  +  Sl^  +  d04^f2  +  S2  ^ 


+  d05f0  +  d06^ri  +  ^1  ^  +  d07^f2  +  Sp )  ) )  Jo  ^01r  ^ 


+  (  ^  i  ^  1  ^  e^d  12^0^1  *  d13^  ^1^2  *  SpSp)  d  14^0^1 


+  d15^rif2  +  S1(l2^^cos0 


I 

i 
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+  +  e^d1210el  +  d13^flsl  ~  eir2^  +  dli4fOsl 


+  dl^(riS2  "  r2S1)))sinO)J1(Snr) 


+  ((d21f2  +  C^d22f0l2  +  d23(fl  Sl)  +  d2^:'of2 


+  d  (f ^  -  g2 ) ) )cos2o 
do  1  1 


+  (d21g1  +  e(d22fOg2  +  2d23flfil  +  d2*T  0S2 


+  2d25f1g1))sin2a)J2(S21r)) 


(5.12) 


The  coefficients  d,,  ,  d,  ,  d^  are  calculated  in  the  same 

Os  Is  2s 

fashion  as  those  discussed  previously  and  the  values  are 
listed  in  table  8  of  the  appendix. 

Some  typical  results  for  wall  pressure  waveforms  are 
presented  in  figures  6  through  37. 

Figures  6  through  9  correspond  to  instantaneous 
combustion  with  e  =  .05,  w  =  .1,  n  =  175,  and  initial 
conditions  (5-10).  This  leads  to  a  stable  standing  wave 
oscillation  and  the  pressure  can  be  observed  to  decrease 
gradually.  However, by  changing  the  interaction  index  to 
n  =  220  (Figures  10-13)  the  pressure  is  seen  to  grow 
gradually.  This  is  an  unstable  situation.  In  both  cases 
the  response  is  dominated  by  the  IT  mode  but  distorted  to 
some  extent  by  the  presence  of  the  2T  and  1R  components. 

Figures  l(t  through  17  correspond  to  history -dependent 


combustion  with  n  =  3^5,  while  the  other  parameters  are  the 
same  as  before.  This  is  a  stable  situation  and  the  amplitude 
of  the  pressure  slowly  decays  with  time. 

Figures  18  through  21  are  obtained  by  changing  the 
interaction  index  to  n  =  352.  This  is  an  unstable  situation 
as  indicated  by  the  growth  of  the  pressure  amplitude  with 
time.  In  both  of  these  cases  the  presence  of  the  2T  and  1R 
components  is  much  more  noticable  than  it  was  in  the  first 
two  situations. 

Figures  22  through  25  correspond  to  instantaneous 
combustion  with  n  =  180,  e  =  .05,  w  =  .1,  and  initial 
condition  (5.11).  This  produces  a  traveling  wave.  The 
pressure  is  observed  to  be  decreasing  with  time  (a  stable 
situation)  while  changing  the  interaction  index  n  to  210 
(Figures  26  through  29)  causes  the  pressure  to  increase 
(an  unstable  situation).  Figures  30  through  37  show  the 
significance  of  the  time  delay  function.  Figures  30  through 
33  are  obtained  by  setting  t  =  it  and  n  =  197.  This  is  a 
stable  case.  The  pressure  is  observed  to  decrease.  Setting 
n  =  199  (Figures  3^  through  37),  on  the  other  hand,  produces 
an  unstable  situation  when  the  pressure  increases  very 
rapidly.  In  both  of  these  situations  the  response  appears 
to  be  dominated  by  2T  contribution  to  the  pressure. 

From  these  figures  one  can  conclude  that  the  pressure 
waveforms  exhibit  a  strong  second  harmonic  distortion  and 
this  distortion  arises  from  the  effect  of  the  quadratic 
nonlinear  terms.  A  variety  of  behaviors  are  possible 


depending  on  the  nature  of  the  combustion  process  and  the 
parametric  values  involved. 
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Chapter  6 

ONE-DIMENSIONAL  MODEL 

In  this  chapter  an  annular  combustion  chamber  with  a 
gap  width  much  smaller  than  the  inner  radius  is  considered. 
This  geometry  will  henceforth  be  referred  to  as  that  of  a 
narrow  annulus.  While  this  geometry  is  not  of  much  practical 
interest, it  is  quite  useful  for  the  purposes  of  analysis. 

This  is  because  only  one  space  coordinate  is  needed  to 
describe  the  problem.  This  makes  a  direct  finite-difference 
numerical  solution  of  the  original  partial  differential 
equation  feasible  and  also  simplifies  the  algebra  required 
to  carry  out  the  Galerkin  modal  analysis.  In  what  follows 
three  questions  will  be  investigated.  First, the  effect  of 
changing  the  numerical  values  of  certain  coefficients 
appearing  in  the  governing  equations  for  the  modal  amplitudes 
will  be  discussed.  Second, the  Galerkin  solution  will  be 
checked  using  a  finite  difference  numerical  solution  of  the 
complete  equation.  Third,  numerical  solutions  of  the 
complete  wave  equation  using  the  vaporization-limited 
burning-rate  law  will  be  compared  to  similar  solutions 
associated  with  the  phenomenological  burning-rate  law 
employed  In  the  previous  chapter. 

The  appropriate  wave  equation  for  transverse  wave 
motion  in  a  narrow  annulus  can  be  obtained  from  (5.6)  by 


setting  r  =  1  and  3r  =  o.  To  further  simplify  the  results 
the  parametic  values  y  =  1  (isothermal  process)  and  J  =  0 
(instantenous  burning  response)  will  be  employed.  For  this 


situation  (5.6)  simplifies  to 

3tt*  +  £at*  ~  D o e ^  +  2*e3e<#>3et^  +  wne(3^)2  =  0.  (6.1) 

'fo  carry  out  the  modal  analysis,  it  is  assumed  that 
the  potential  function  can  be  expressed  as 

*  "  ^(tJcosG  +  f^  (t  )  cos  20  +  g1(t)sino  +g.,(t  )sln2e.  (6.2) 


where  r.^  =  l,  =  2,  A^  =  2,  A-,  =  2,  A.  =  1  and  A^  =  .5. 

The  symbols  il0,  A^,  A.,,  A.,  and  A^  have  been  Inserted  to 

illustrate  the  effect  of  changing  their  numerical  values. 

McDonald  (7)  in  his  analysis  of  combustion- 
instability  in  an  annulus  found  that  for  a  given  value  of  e 
the  value  of  n  required  to  produce  instability  was  approx¬ 
imately  twice  as  high  for  standing  waves  as  for  traveling 
waves.  The  results  discussed  in  the  previous  chapter  for 
a  full  cylinder  indicated  no  such  relationship.  Equations 
(v.2)  were  employed  in  an  attempt  to  determine  the  factors 
which  have  a  significant  effect  on  the  stability  boundary. 

Several  calculations  were  made  and  a  representative  sample 
of  the  data  thus  obtained  is  presented  in  Tables  1  and  2. 

It  was  determined  by  McDonald  that  the  terms  representing 
gas-dynamic  nonlinearities,  had  a  small  qualitative  effect 
on  stability  calculations.  Thus  the  coefficients  An  and  A_, 

1  J) 

were  held  fixed  during  these  calculations. 

The  entries  in  the  first  two  lines  of  each  table 
were  computed  using  the  correct  equation  for  an  annulus.  It 
can  be  seen  that  the  standing  wave  is  twice  as  stable  as  the 
traveling  wave,  in  agreement  with  the  results  of  McDonald. 

The  third  and  fourth  lines  in  each  table  were  computed  by 
changing  A;.  from  .5  to  .77  (This  makes  the  ratio  A^/A 
the  same  as  the  ratio  of  the  corresponding  terms  in  the 
governing  equations  for  the  full  cylinder.).  It  can  be 
seen  that  this  lowers  the  stability  limit  in  all  cases  but 

*  . 
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does  not  alter  the  fact  that  an  initial  disturbance  in  the 
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form  of  a  standing  wave  is  twice  as  stable  as  one  in  the 
fora  of  a  traveling  wave. 

The  fifth  and  sixth  lines  in  Tables  1  and  2  are  found 
by  restoring  ;V  to  its  original  value  .5  and  changing 
from  2.00  to  1.66  (This  rakes  the  ratio  0,  equal  to  the 

t_  X 

ratio  S-,,/2.^  associated  with  the  full  cylinder.).  For  this 
situation  it  can  be  seen  that  the  standing  wave  is  approx¬ 
imately  ten  times  as  stable  as  the  traveling  wave.  Further¬ 
more,  the  effect  of  the  value  of  sly  on  the  l°oatior,  of 
stability  boundary  is  much  greater  for  a  standing-wave 
motion  than  for  a  traveling-wave  motion.  The  last  two  lines 
are  associated  with  implementing  both  of  the  changes  dis¬ 
cussed  ibove  simultaneously.  These  results  confirm  that  the 
change  in  A  i(  lowers  the  stability  limit  in  all  cases  but 
toes  not  affect  the  relative  stability  of  standing-wave  and 
t raveling-wave  disturbances . 

The  data  presented  above  indicate  that  standing  waves 
will  be  twice  as  stable  as  traveling  waves  only  under  very 
special  circumstances  ( n 0 / U ^  *  2).  There  is  no  reason  to 
expe'i  this  to  be  a  characteristic  of  other  systems  exhibit¬ 
ing  combustion  instability  as,  in  fact,  it  is  not  for  a  full 
cylinder . 

In  this  section,  a  finite  different  procedure  is 
employed  to  generate  a  set  of  second  order  differential 
equations.  The  central  difference  formulas 


V,  ■  <"L+1  - 
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3ee»i  =  (»i  +  i  '  2*i  +H-i)/(A0)2  (6>U) 

are  used  to  produce  the  following  set  of  governing  equations 
$.  =  — w>>,  +  +  ^  -  2<p^  +  _  c  (  y-1 )  ^^  )/ (  AS  )  2 

-  c(*i+i  -  ‘  *i-i,/(2ie)2 

_  eWi  2  <  i  <  N  -  1  (6.5) 

where  N  is  the  number  of  points  employed.  A  forward 
difference  formula 

8  =  (-3*,  +  4<}>  -♦  )/(2A0)  (6.6) 

el  l  2  a 

and  backward  difference  formula 


30<f>N  "  ^N-2  "  44>N-1  +  34>n)/(2&0) 


(6.7) 


are  used  for  the  boundary  conditions 


30*(O)  =  0,  3e4>(")  =  0. 

Hence,  along  the  boundaries,  (6.5)  becomes 

=  -w$_  +  Z*  ( 4>  0  -  4>.,)(1  -  c(  Y~l)  v’0)/(3ao2) 


(6.8) 


t 


-  8  G  C  <J>  _  -  <?„)(<>_  -<?>,  )/(9A02  )  -  £  w 

J  <-  J  tv  c. 

and  =  -w^N_1  +  2($n_2  "  )  ( 1-c ( y-D ❖VJ_1 ) /(  3A0 2  ) 

-  8e(«N-1  -  -  *m.2)/(948Z>  -  (6'9> 

For  the  phenomenological  model,  the  burning  rate 
function  becomes 

w,  =  nw($1+1  -  4>1_1)2/C4&92)  2  <  i  <  N- 1 .  (6.10) 

Combining  (6.8)  and  (6.10)  one  obtains 


w.  =  w((l  -e<j>i/2)(l  +  (<!>1  +  1 

-  ♦.  )2/(2A9ut )2)k  -  1 )/ e2  2  <  i  <  N-l  (6.12) 


and, along  the  boundary,  yields 
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w.  =  w((l  -  JjeJ_)(l  +  (2  C  <?  _  -  a  )/(3A8u,  ) )  2  )’6  -  l)/£2 

c  J  C  J-J 

ww_l=w((1-HeiN_l)(1+C2(«N_1-<pN_2)/(3A6uL))2)i4-l)/E2,  (6.13) 

accordingly . 

Using  either  the  phenomono logical  or  the  vaporization- 
limited  combustion  model, (6. 5)  and  (6.9)  can  be  solved  by 
the  Runge-Kutta  method  to  obtain  $ . ,  i  =  2,3,-..,  N-l. 

Then  the  modal  amplitudes  can  be  obtained  by  using  the 
Fourier-cosine  transform 


7T 

f.(t)  =2/  $cosi0de/n.  (6.14) 

i  o 

These  integrals  must  be  computed  numerically  since  <>  is 
known  only  at  the  grid  points. 

Some  typical  results  for  standing  waves  are  presented 
in  Table  3-  They  are  intended  to  illustrate  the 
accuracy  of  the  two-term  Galerkin  solution  and  to  compare 
the  results  associated  with  the  phenomenological  and 
vaporization-limited  combustion  laws.  The  column  labeled 
PI  indicates  the  results  with  a  two-term  Galerkin  solution 
using  the  phenomenological  model,  the  column  labeled  PF 
contains  the  results  of  a  finite- difference  solution  of 
these  equations,  and  the  column  labeled  VF  presents  results 
of  a  finite  difference  solution  using  the  vaporization- 


limited  model. 
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Table  3 


Stability  Limits  For  Phenomenological  and 
Vaporization-Limited  Combustion  Model 


j= - 

PG 

PF 

VF 

e 

n 

n 

UL 

•  1 

46 

28 

1.3 

.05 

91 

56 

2 . 1 

.01 

_ i 

45  2 

286 

5.2 

It  can  be  seen  that  the  order  of  magnitude  of  the 
interaction  index  n  required  for  instability  for  both  the 
finite  difference  and  Galerkin  methods  are  roughly  the  same 
but  the  stability  boundaries  predicted  by  the  Galerkin  method 
are  approximately  twice  as  high  as  those  predicted  by  the 
Finite-difference  method.  It  is  to  be  expected  that  the 
Galerkin  method  will  lead  to  higher  stability  limits  than 
the  use  of  an  exact  solution  procedure.  This  can  be  explained 
as  follows.  The  instability  mechanism  is  basically  a  feed¬ 
back  process.  Consider  the  form  of  (6.3)  associated  with 
standing  v/aves  in  an  annulus .  This  is 

fl  +  -fl  +  fl  +  2e(f1f‘2  +  f2fl)  +  2£wnf1f2  =  0 

f2  +  wf2  +  4f  -  ef1f1  -  *sewnfj  =  0.  (6.14) 


non  1 1  no an ! t  1  os 


For  clarity  the  gas -dynamic 
to  give 


will  be  neglected 


l'j  +  wfj  +  f  j  +  2  fKn f  2 1' o  =  0 


f +  wfn  +  -’I  f 2  -  Sewnf' 


(6.  rd 


A  one-term  tialerkln  solution  would  lead  to  the 


e  iuat  ion 


r ,  +  w  f  -j  +  r  =  o . 


This  would  indicate  uncond it ional  stability.  Instability 
can  arise  only  when  the  presence  of  the  last  term  in  (6.15b'1 
causes  f,  to  grow  and  this  then  causes  f^  to  grow  because 
of  the  last  term  in  (6.15a).  The  growth  of  t'_,  also  produces 
the  growth  of  higher  modes  which  in  turn  causes  additional 
growth  of  fj.  The  contrlbut ions  of  these  higher  modes  are 


neglected  in  a  two-term  ualerkln  analysis  and  thus  the 


energy 


Input  due  to  unsteady  burning  is  underestimated.  For  this 
reason  the  d.ilerktn  analysis  will  overestimate  the  stability 
of  the  system.  This  overestimat ion  will  decrease  as  the 
number  of  terms  retained  is  increased.  An  example  of  this 
can  be  seen  by  comparing  the  one-  and  two-term  analyses. 

A  one-term  solution  predicts,  tint  the  system  Is  always 
.-.•able.  A  two-term  solution  predicts  the  correct  gualita- 
t  i  ve  behavior  but  overest  lnat.es  the  system's  stabl  ltty  by 


roughly  a  factor  of  two.  It  is  to  be  expected  that  further 
increases  in  accuracy  can  be  obtained  by  keeping  more  terms 
in  the  assumed  solution  but  that  this  will  not  seriously 
affect  the  quatitative  prediction. 

In  Chapter  5,  it  was  found  that  equating  dw/dur2  the 

U 

phenomenological  and  vaporization-limited  combustion  laws 
at  u*  =  0  lead  to  the  formula 

n  =  l/(4e2u2)  .  (6.16) 

Assuming  that  the  actual  data  can  be  fitted  to  a  formula  of 
the  form 

n  =  C/ ( 4  e 2u 2 )  (6.17) 

the  data  give  in  Table  3  provide  an  opportunity  to  estimate 
C.  The  results  are  shown  below. 


_ L_ 

.  I 

.05 

.01 

_ c_ 

*—• 

CO 

o 

2 . 47 

3.09  1 

It  appears  that  a  value  of  C  =  2.5  would  give  accept¬ 
able  accuracy  over  this  range  of  e. 

For  a  given  value  of  ur  the  n  computed  from  (6.16) 
will  overestimate  the  burning  rate.  Thus  it  might  be 
thought  that  C  should  be  less  than  unity.  The  stability 
criterion  was,  however,  based  on  the  magnitude  of  the  modal 
amplitudes.  Inspection  of  Tables  4  and  5  shows  that  the 


f 


vaporization-limited  combustion  model  Involves  the  higher 
modes  to  a  greater  extent  than  does  the  phenomenological 
combustion  model.  Thus, a  given  value  of  Uj  can  be  associated 
with  lower  Individual  modal  amplitudes  in  the  former  case 
than  in  the  latter.  These  two  effects  Interact  and  calci la- 
tion  shows  that  C  is  greater  than  unity  in  this  range  of  e. 
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TA8LE  A 

MODAL  AMPLI  TUOE.  FOR  PHF  NO  M£  ML  CG  I  C  AL  MODEL 
IN  ANNLLAR  C0M8GST0R 


T 

RODE  1 

MODE  2 

MODE  3 

MODE  4 

NODE  5 

MODE  6 

MODE  7 

MODE  3 

0.  CO 

1  .  OCu 

0 . 000 

-0.000 

0.000 

-0.000 

0.000 

-0. 000 

0.  CO  0 

1.25 

0.  3  37 

0.  C  2  9 

C.  002 

0 .000 

0.000 

-0.000 

-0 . 000 

-0.000 

2.50 

*0.657 

-0.018 

-0. 000 

0.000 

0.  GCO 

-C.003 

0.  000 

O.COO 

3.  75 

”0- 7C0 

0 .064 

-0.003 

0.001 

-0.00  0 

0.  ooc 

0.  000 

-0.000 

5.  CO 

0.  2 C  9 

-0.  04  3 

O.OC9 

-0.001 

-o  .OCO 

G  .000 

-0.000 

0.  ooc 

6  •  c  5 

0.751 

0.040 

-0.006 

-0.003 

-o. 0C1 

0.00  0 

0.000 

0-000 

7.  50 

0.237 

0.025 

0.  OCO 

-0.002 

-0.0C1 

- C. 00  0 

-0.000 

-0. 000 

3.75 

*0.535 

-0 . C60 

C.015 

0.004 

-0.002 

-0.000 

0 . 000 

0.000 

1  C.  C  0 

*0.523 

0.105 

-0.023 

0.009 

-O.OC  3 

0.001 

-3. 000 

0.000 

11.25 

0.137 

-0.038 

0.  024 

-0.000 

-0. CC  3 

0.302 

-0.000 

-0.000 

12.50 

0.555 

0  .  C62 

-0.009 

-0.011 

-0.C05 

-0.001 

0.001 

0.  001 

1  3.  75 

0.153 

0.011 

-0.010 

-0.  310 

-0.006 

-0.003 

-3.002 

-0.00  1 

1  5.  CO 

-0. 4A0 

-0 .065 

0.0  32 

0.007 

-o. 0C7 

-0.301 

0.002 

o.ccc 

1  6  .  c  5 

-0.390 

0.112 

-C.043 

0.020 

-0.010 

0.  00  5 

-0. 002 

0.  00  1 

17.50 

0  .  1  S  3 

-0 . 105 

0.0  53 

0.00  5 

-0  .010 

0  . 0  G  4 

3  .  C  0  1 

-0.002 

15.75 

0  .  <.  c  3 

0.073 

-  0 . 004 

-  0  .  C  1  7 

-0.013 

-0.004 

0.031 

0.  002 

2C.C0 

0.055 

-0. 009 

-0. 023 

-0.020 

-0.014 

-0-008 

-0. 004 

-0. 002 

2  1.25 

-0.  3 7  A 

-0 .050 

0.0  45 

0.034 

-0.015 

u  .00  I 

3.006 

-0.000 

2  2 . 5  0 

-0.2  J  7 

0.101 

-0.043 

0.02  5 

-0.013 

0.011 

-0.  007 

0.033 

2  3.  7  5 

0.  1  5  A 

-0.107 

0.025 

0.015 

-0.013 

0.  005 

0-004 

-0.005 

25.  CO 

0.37  1 

0.089 

0.009 

“0.01b 

-0.019 

-0 .c 10 

0.00! 

0.  00  5 

2  6.25 

0.0  50 

-0.  C  31 

-  0.  0  37 

-  0.  030 

-0  .020 

-0.012 

-0.006 

-0.002 

2  7.50 

-0.  5  2  7 

-0.027 

0.05? 

-0.005 

-0.021 

C.0C5 

0.010 

-0- 002 

Jf>.  75 

-0. 2Cb 

0.031 

-0.0 A4 

0  .0  30 

-0.023 

0.  017 

-0.  oil 

0.  006 

5  0.  CO 

0.191 

-0.102 

0.013 

0.023 

-0.02  5 

0.003 

0 .008 

-0.006 

3  1.25 

0  •  3C2 

0 . 098 

0.026 

-0.009 

-0.C23 

-0.01/ 

-0.003 

0  •  0  0  h 

3  2.  50 

-0.019 

-0 . 054 

-0.  051 

-0. C33 

-0. 024 

-0.012 

-0.004 

0.00  1 

3  3.75 

-0.297 

0  .  coc 

C.C56 

-0.019 

-0.024 

0.0  12 

O.C  10 

-0.005 

35.  CO 

-0.  1  A  1 

0.053 

-0.034 

0.027 

-  0. 02  A 

0.  020 

-0.013 

0.008 

36.  25 

0.  2  07 

-C . 093 

-0.005 

0.040 

-0.023 

-0.003 

0.012 

-0.005 

3  7.50 

0.  252 

0.105 

0.046 

0.004 

-0.023 

-0  .322 

-0.008 

0.  001 

3  6.75 

-0.067 

-0. 079 

-0.064 

-0.  0A2 

-0.02? 

-0  .007 

0.002 

0.005 

<.0.00 

-0.2  5  A 

0.031 

0.056 

-0.03b 

-0.020 

0.019 

0.005 

-0.00  8 

<.1.25 

-  C  .  0  6  5 

0.  C33 

-0.019 

0.018 

-0.019 

0.016 

-0.010 

0-  00  7 

<•2.50 

C.  236 

-0.087 

-0.029 

0.052 

-0.017 

-0.013 

0.015 

-0.001 

<.3.75 

C  .2  18 

0.120 

0.070 

0.023 

-0. 01  A 

-0-022 

-0.014 

-0. OCfa 

<.5.  CO 

-0.12? 

-0.109 

-0.076 

-0.039 

-0-011 

0.006 

0.012 

0.010 

A  6  .  2  5 

-0.292 

0  .  C  70 

0.0  A9 

-0.055 

-0.0 07 

0.025 

-0.007 

-0.008 

A  7.50 

-0 .02  9 

-0.001 

C.005 

-0.001 

-O.OC3 

0.  002 

0.  000 

-0.000 

Ae.  75 

0.  290 

-0.073 

-0.063 

0.058 

0.  OC  3 

-  0.  02b 

0.007 

0.  0 u  9 

50.  CO 

0.192 

0.135 

0.0  98 

0.053 

0.0  11 

-0.0  0  8 

-0.0  1  A 

-c.  o;  3 

5  1.25 

-0.207 

-0.152 

-  0.  0  81 

-0. 020 

0.016 

0.026 

0.019 

0.006 

52.50 

-0-328 

0.134 

0.022 

-0.07  3 

0.025 

0.015 

-0.02? 

0.00  6 

5  3.  75 

0.055 

-0.060 

0.0  55 

-0 .0 A4 

0 .0  3  A 

-  0. 02e 

0.  02? 

-0.015 

55.  CO 

0. 4  C  A 

-0 . C  33 

-0. 1  19 

0.037 

O.OAb 

-0.022 

-0.020 

0.  Cl  2 

5  6.25 

0.151 

0.140 

0.125 

0.094 

0.  061 

C.  0  38 

0.019 

0.008 

5  7.50 

-0.AC3 

-0.219 

-0.0 A5 

0.054 

0.065 

0.025 

-0. 01 3 

-0.  022 

5  6.75 

-0  .  A  1  2 

0.  ?7o 

-0.092 

-0.036 

0.071 

-0.05  1 

3.013 

0.01  2 

6C-  CO 

0.324 

-0.246 

0.193 

-0.117 

0.060 

-0.017 

-0.011 

0.  022 
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TABLE  5 


MODAL  amplitude  fop  vaporization  model 

IN  ANNULAR  COMOUSTCR 


T 

C  JOE  l 

MODE  2 

MODE  3 

MODE  4 

MODE  5 

M  UD  £  6 

MODE  7 

MODE  8 

c.co 

1 . 0  C  0 

0  .  OCO 

-o.coo 

0.000 

-0 .006 

0  .0  00 

-0. 000 

0.  COO 

1  .28 

0.  394 

-0.  1  93 

0.  0  39 

-0.007 

0.002 

- 0  .0  0  1 

0  .  OilO 

-0.000 

2.50 

*0.455 

0 .095 

-0.033 

0.009 

-0.  0C2 

0 . 0  0  u 

0.  000 

-0.000 

5.  75 

-0.514 

-0.111 

0.042 

O.C  10 

-0.00 7 

-  0. 00  3 

0.  002 

0.  000 

5.  CC 

0.  06  0 

-C.  C35 

-0.0  36 

-0  .029 

-0  .020 

-0  .0  12 

-0.007 

-0. CO  3 

6.25 

0.443 

0.106 

-0.022 

-0.041 

”0.021 

-  0.002 

0.006 

0.00  6 

7.  50 

0.  25  5 

-0.225 

0. 077 

0.020 

-0.043 

0.019 

0.007 

-0.01  5 

d.  75 

-0.  524 

0.160 

-0.104 

0  .  C  7  1 

-0.049 

0.034 

-0.024 

0.017 

1  C.CO 

-0.717 

-0.029 

0.192 

-0.037 

-0. G72 

0.04  ? 

0.020 

-0.032 

11.25 

-0.224 

-0.455 

-0. 330 

-0.172 

-0. 040 

0.34  1 

0.068 

0.04  7 

1  2.50 

1.291 

0. 558 

C.  300 

0.131 

0.040 

-0.003 

-0.021 

-0.  030 

13.75 

2.175 

-0.505 

-0.474 

0.  00  7 

0.172 

0.024 

-0-070 

-0.05S 

1  5. CO 

0 . 2  o  4 

-0.967 

0.644 

-0.327 

C.  068 

0-07  6 

-0.105 

0.C5  7 

16.28 

-2.  7  C  7 

0.  857 

-G.  311 

0.100 

-0.053 

0.00  5 

0. 01  9 

-0. C  5 1 

17.60 

-5.416 

-1.125 

0.  362 

O.C  37 

-0.195 

0  .0  5  2 

0.072 

-0. 0  7  4 

18.75 

0  •  4  4  d 

-0 .794 

-0.551 

-0.324 

-0.057 

0 . 0  C  5 

0.110 

0.  C56 

2  C.  CO 

3.619 

0  .  54  o 

0.  200 

0. 07  7 

0.GC1 

-  0.  024 

-0.039 

-0.044 

2  1.25 

2.747 

-1.259 

-0.224 

0.106 

0  .  l6o 

-C  .02  4 

-0.087 

-0.025 

2  2.50 

*1.442 

-C . 323 

0.  340 

-0-27  1 

0.  145 

-0.042 

-0.  0  52 

0. 06  3 

2  3.  7  5 

-3.  22  9 

0.136 

-0.115 

“0.060 

0.  063 

-0.030 

0.058 

-0-0.58 

2  8  .  C  0 

-1.  52  9 

-1-066 

0.005 

0 . 10  1 

-0.028 

-0.10  2 

0.0  12 

0.  05q 

2  6.28 

1  .  794 

0.  C  24 

-  0.  009 

-  0.  005 

-0.08a 

-0.083 

-0.067 

-0.049 

2  7.50 

2.395 

-0.063 

-0-093 

-0.167 

-0- C39 

-C.030 

0-022 

0.04  1 

2  6.75 

0  .  4  t  3 

-0.831 

0.2  35 

0  .090 

-0.11  4 

-  0.  006 

0.  0  69 

-  0. 02  7 

3  C.  CG 

-1  .  3C7 

0.272 

-0.198 

0.095 

-0.067 

0.045 

-0.031 

0.021 

31.25 

-2.093 

-0 .294 

0.261 

-0.14  5 

-0. C24 

0.07  5 

-0.082 

-0. 00 3 

3  9.  60 

-0.,  06  5 

-0.796 

-  0. 404 

-0.100 

C.  077 

0.39  3 

0.025 

-0.040 

3  3.75 

2.115 

0.4  17 

0.2  50 

0.120 

0.079 

0.049 

0.033 

0.021 

35.  CO 

2.172 

-0 .647 

-0. 302 

-0.013 

0.123 

0.062 

-0.032 

-0.059 

3  6.  25 

-0.  447 

-0.603 

0.422 

-0.241 

0-058 

0.049 

-0.079 

0.  050 

3  7.50 

-2. 596 

0. 434 

-0.193 

0.046 

-0  .004 

-0.025 

0.035 

-0.  0  39 

3  °.  75 

-2.044 

-1 . 000 

0.179 

0-126 

-0.127 

-0.044 

0.074 

O.OOC 

MO. CO 

1.142 

-0.253 

-0.260 

-0.241 

-0.154 

-C.072 

-0 .004 

0.038 

A  1  .  2  5 

2.34  3 

0.253 

0.0  55 

-0  .083 

-0  .089 

-0.08  1 

-0.  054 

-0.02  7 

A  2.  50 

1.4  97 

-1.140 

0.032 

0.196 

0  .007 

-0.109 

0.006 

0  .  C  6  5 

4  5.75 

*  1  -  7  2  1 

0.032 

0.011 

-0.091 

0.  095 

-0.09  3 

0.077 

-0. 057 

4  5.  CC 

-2.759 

-0.014 

0.116 

-0.183 

0.105 

-  0. 0  38 

-0. 020 

0.  04  4 

4  5.25 

-0.007 

-1.004 

-0.251 

0.12b 

0.118 

-0 .0  3  3 

-0.00  1 

-0.012 

47.50 

2.050 

0.301 

0.193 

0.08  3 

0.050 

0.024 

0.  010 

-0. 001 

4  9.75 

2.  51  4 

-0.289 

-0-260 

-0.186 

-0. Cll 

0.  062 

0.  068 

0.02  3 

5  0.  CO 

0 . 2  r'  b 

-0.932 

0.405 

-0.043 

-0.117 

0.082 

0.014 

-0.  Oo  3 

6  1.25 

-2.255 

0.  39  7 

-0.  263 

0.134 

-0.096 

0.066 

-0.051 

0.04  0 

52.50 

*2.237 

-0.575 

0.  312 

-0. 076 

-0.106 

0.034 

-0. 002 

-0. 05  3 

5  3.  7  j 

0.3  31 

-0 . 709 

-0.439 

-0.200 

-0.00  0 

0.08  3 

0.  07  4 

C.  01  7 

55.  CC 

2.477 

0.407 

0.229 

0.032 

0.038 

0  .00  4 

-0.011 

-0.021 

5  6.25 

2 .025 

-0 . 053 

-0.236 

0.07  3 

0.141 

-  0.  002 

-0.0/2 

-0.0  31 

5  7.  c0 

-0.9  15 

-0.377 

0.  322 

-0.252 

0.129 

-0.03C 

-0.0  36 

0.  060 

5  9.  75 

-2. 669 

0  .  3  0  3 

-0.109 

-0.033 

0  .057 

-  0  •  0  b  8 

0.057 

-0.042 

6C.  CO 

-1.593 

-1.060 

0.033 

0.179 

-0.051 

-0.  099 

0.028 

0.055 

-  >'  unn.imw i 


A. 


Chapter  7 


CONCLUSIONS 

The  primary  objective  of  this  study  was  the 
development  of  a  new  analytical  technique  to  be  used  in  the 
solution  of  nonlinear  velocity-sensitive  combustion 
instability  problems.  Such  a  method  should  be  relatively 
easy  to  apply  and  should  require  relatively  little 
computation  time. 

In  an  attempt  to  achieve  this  aim,  the  orthogonal 
collocation  method  was  investigated  first.  However,  it  was 
found  that  the  results  were  heavily  dependent  on  the  location 
of  the  collocation  points  and  characteristics  of  the 
equations.  Therefore,  the  method  was  rejected  as  unreliable. 

Next,  the  Galerkin  method,  which  has  proved  to  be 
very  successful  in  analysis  of  the  pressure  sensitive 
combustion  Instability,  was  considered.  This  method  proved 
to  work  very  well.  It  was  found  that  the  pressure  waveforms 
exhibit  a  strong  second  harmonic  distortion  and  a  variety 
of  behaviors  are  possible  depending  on  the  nature  of  the 
combustion  process  and  the  parametric  values  involved. 

Finally,  a  one-dimensional  model  provided  further 
insight  into  the  problem  by  allowing  a  comparison  of 
Galerkin  solutions  with  more  exact  finite-difference 
computations . 
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specific  assumptions  used  to  derive  it.  This  is  demonstrated 
by  the  fact  that  Powell  (1)  developed  an  equation  of  a  very 
similar  form  using  a  different  set  of  assumptions.  (2)  The 
orthogonal  collocation  method  is  unsuited  to  solution  of 
problems  of  the  type  under  discussion  here.  (3)  The  Galerkin 
method  is  well  suited  to  the  solution  of  such  problems.  (*0 
ftability  boundaries  and  pressure  wave  forms  appear  to  be 
highly  dependent  on  the  parameters  of  the  problem  (interac¬ 
tion  index,  time  delay,  steady-state  burning  rate,  etc.). 
Furthermore , there  appears  to  be  no  clear  pattern  to  the 
computed  results.  (5)  The  phenomenological  burning-rate  law 
employed  in  the  majority  of  the  work  discussed  predicts 
results  which  are  qualitatively  similar  to  those  associated 
with  the  vaporisation-limited  burning-rate  law  used  by 
previous  investigators.  The  phenomenological  law,  further¬ 
more,  can  be  used  in  conjunction  with  the  Galerkin  method 
while  the  vaporization-limited  law  cannot.  (6)  The  computer 
programs  developed  in  the  course  of  this  work  can  be 
employed  to  determine  stability  boundaries  and  pressure 
waveforms  without  the  expenditure  of  excessive  computer 
time.  This  is  important  because  the  lack  of  clear  trends 
discussed  above  make  an  individual  analysis  of  each 
situation  desirable. 
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The  Wall  Pressure  Waveforms  for  Traveling  Waves 
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APPENDIX  A 


DERIVATION  OF  THE  COEFFICIENTS  \ 

In  this  section,  the  coefficients  b,  ,  ,  C  ,  and  d, 

v)  "*■  J  O 

appearing  in  chapter  3  and  5  are  derived  and  their  numerical 
values  are  presented  in  Table  6,  7, and  8  respectively. 

In  chapter  3,  for  initial  condition  (3.7),  one  assumes 


3, ,.$-5  -  v2$0  =  T.  H,  (r)cos(me)sin(231  _t) . 

^  L  tL  <-  HI—  0  m  -A- 

By  comparing  coefficient  yields 


H0(r>  - 


v-  (  c  2  t2 


-  ^5. 


reJi/r 


4J0/r: 


H1(r)  =  0 

H2(r)  =  ^Sn(2S21J2  -  (y-DS^jJ  -  4S11JQJ1/r) 

H  (r)  =  0  for  m=  3,  4,...,«>. 

Then  expanding  in  a  Fourier-Bessel  series  and  expanding 

H  (r)  in  a  Bessel  series  leads  to 
m 
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This  integral  must  be  computed  numerically.  For  the  sake 
of  brevity  only  the  first  five  coefficients  in  each  series 
are  calculated  and  presented  in  Table  6. 


Table  6 


5ESSEL  SERIES  COEFFICIENTS  FOR  ANALYICAL  SOLUTION 


n 

b0n 

b2n 

0 

1. 11738-0. 37246 y 

1 

0.50462+0. 37912y 

0. 27181-1. 09482y 

2 

-0. 09246-0. 00802y 

-0. 11390-0. OlllOy 

3 

0. 05019+0. 00134y 

0. 05479+0. 00208y 

A 

-0. 03290-0. 00067y 

-0. 03451-0. 00071y 

5 

i 

I 

0 . 02  375+0 . 00030 y 

L  -  - — 

0. 02451+0. 00033y 

1 

In  Chapter  5,  equation  (5*8)  can  be  written  as 


*  =  Vo  +  Vl  +  Vs  +  sl*3  +  *2*4 


( A—  1 ) 


where 


=  J, 


<*>  i  =  JjCosO,  <t>  =  J0cos26 


=  J^sine,  =  J^sinSe. 


By  substituting  (A-l)  into  the  governing  equation  (5.6), 
yields 


R  =  D(*) 


(  A-2 ) 
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where  D  is  the  nonlinear  differential  operator  of  equation 
(5.6).  If  ( A-l)  represented  the  exact  solution  to  (5.6),  R 
would  vanish.  Since  this  is  not  the  case  R  has  a  finite 
value.  The  Galerkin  procedure  consists  of  making  R  orthog¬ 
onal  to  each  of  the  <t^'s.  This  leads  to  the  equations 

1  2  IT 

/  /  R$  rdrd  =0  i  =  0,1,. ..,4  (A-3) 

0  0  i 

and  a  set  of  five  second  order  differential  equations  are 
generated.  The  angular  part  of  integration  can  be  performed 
In  closed  form  while  the  radial  part  of  integration  must  be 
computed  numerically  and  leads  to  the  integral 

=  6k/QPi(r) Jk(Sklr)rdr  (A-4) 

™ere  Bk  =  SS^/iCSj^  -  k2 )  CSkl)  } .  (A-5) 

The  numerical  values  of  are  listed  in  Table  7  for 
Y  =  1.2.  Also  listed  in  this  table  are  the  functional  forms 
of  F, 's.  In  writing  these  the  notation 

Ri  =  3ii(Y-1)Ji  +  SilJi/r  "  l2jj/r2>  1=0, 1,2  ( A-6 ) 

is  used.  Similarly  the  numerical  values  of 

dij  =  6i/J0GijJiCSilr)rdr  (A-7) 
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Table  7 


Table  8 


COEFFICIENTS  APPEARING  IN  EQUATION  (5.12) 


i  1 

1.000 

Ch 

O 

Tj> 

o 

2  1 
! 

1.29  3 

1  c- ?  y  l  2 
•^0  0 

3  ; 

0.21)0 

^(S^jJ'j2  + 

Jf/r2) 

r 

-0.09S 

H(S2jJJ2/2 

+  2J2/r2) 

5 

-0.176 

7  2 

0 

6 

0.061 

-Vf 

7  j 

0 .049 

_V  J  2 

2 

1 

1.000 

Ji(Snr) 

2  1 

-1.212 

S0 1S 1 1 J o J 1 

0.«7 

h (S !  1S21J  [J 

\  +  2JxJ2/r2) 

4 

0.165 

-J  0J  J 

5 

-0.199 

— I5J  jJ  2 

i 

1.000 

J 2 Cs 2  l^* ) 

2 

'  -1.741 

S01S2lJJJ2 

rjnfinnnooonnnnoonnnoronnnf 
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iSET  L  i  i  T  FREE 

KIN-WING  WONG  JUS  30  15/g 

COMBUSTION  l  NS  T  ABILITY 
w  E I H  QC  GALE  RK I N  MET  HOO 

THIS  PROGRAM  GENERATE  THE  STABILITY  CURVE. 

INPUT  INFORMATION 

T  RE  PROGRAM  USEC  THE  FORMAT  FREE  INPUT  OUTPUT •  A 
CCMMA  SERVES  AS  A  DELIMITER  ANC  the  INPUT  VARIABLES 
ARE  ENTER  IN  THE  FOLLOWING  ORDER 
II  I Nl TI AL  TI ME 

TF  FINAL  TIME 

OT  STEP  SIZE  FOR  TIME 

BURNING  RATE 
TDELAY  DELAY  TIME 

El  I  N  I  TI  AL  EP  S 

EF  FINAL  EPS 

OE  STEP  SIZE  FOR  EPS 

DA  ASSUME  POSI TION  FOR  N 

F  0 » F  0  ARE  THE  INITIAL  CONDITION  MODIFIER  AND 

HAVING  THE  VALUE  EITHER  1  OR  C 


C  ANOTHER  SETS  OF  CATA  CAN  BE  ENTER  3Y  ENTERING 

C  £I»EF»OE»DA,FO  ANO  FD  . 

C 

C I  ME  NS  1  ON  DY(10)»F3(5)»DF(5)»G(10»500  )»Y(  10) 

CCMMQN  EP,  RK , *0 ,EP ANV 
COMMON  /SLIST/S0#S1»S2 

COMMON  ✓BLK1/C1*C2»C3»CA»C5»CS»C7»CS»C9»C18»C11»C12 
1*0 3  *  C 1  4  *  C15*Cl£*Ci7 

LOGICAL  PASS *G ROW 
CATA  MAXITR,EPS/2C*.l/ 

CATA  FI/3. 1A  lb 9/ 

REA0C5*/)  T1*TF*DT.WB, TDELAY 
FRIN  T*  /  /*  T  I  *  TF  *  0  T  »  **8*  TDELAY 
TDELAY  =P I  *  T  CEL  A  Y 
K=TDELAY/DT 
N  1  =  K  ♦  1 

IF  (Kl.GT.5C0)  GO  TO  105 
IF  (K.LT.l.AND.NCF.EC.l  )  K=1 
NSTEP-(TF~TI  )/D T 
M  10 

HHALF  =  0  T  *.5 
N FRO  3  =  0 

9  REA  DC  5*/*£ND=9S9  )  E  I  *  E  F  *  Q£  *  DA  *  F  0  »  CF 
N  F  R  0  9  =  NPRC  8 ♦  1 


4 
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F  F  I  NT  •  //  »N  PF  08 
FUfiT  *//*F0.DF 
NSTEPMEF-EI  >/OE 
EF=E  I 

IF  (EP.LT.O)  E  P  =  .0  1 
AM  =  DA 

FASS=. FALSE. 

CC  IOC  1STEP  =0*NSJEP 
CC  80  ITR=1.MAXITR 
E  F  AN  V  =  EP  *  AN  V 

IF  (EPANV.EC.O  £P ANV  = ANV 
CC  10  1=1. K 
CC  10  j  =  1.N 
1C  G  (  J  .  I  )  =  0  . 

CC  20  1  =  1. N 
2C  1 (I  )=C. 

C 

C  initial  condition 

c 

y<2)=fo(1)Jy<(.>  =  fc(2);y(6>=fq(3);y(6>--fo<<.>;y(10)=fo 

1(5) 

Y (1 >  =  CF(1 )*S0>  Y ( 2) =DF( 2  )  *Sl; Y (5)  =CF  C  3  )*S2 

Y(7)=CF(4)*S1;y<9)=0F(5)*S2 

CALL  CLYFUN  f N . Y » G . A , EPS . K 1 > 

I  =  T  I 
C 

CC  40  I=1.NSTEP 

CALL  RKDELYCN.T  »Y .DY.DT »  G .K  *E  P.4110.HHALF  »Kl  ) 

CC  30  J=Z.N.Z 

1FCA3S(Y(J  )).GT.2)  GO  TO  60 

3  C  CCNTINLE 

4  C  CCNTINUE 

GF0W=  -FALSE  . 

IF  (PASS)  GO  TO  50 
A  N  VD  WN  =  AN  V 
AAV=AN V*DA 
GC  TO  80 

5  G  CCNTINUE 

ANVOWN=ANV 
GC  TO  70 
60  CCNTINUE 

PASS=. TRUE. 

GFOW  =  .TRUE  . 

ANVUP=  ANV 

7  C  CCNTINLE 

IF  ( AeS(ANVLP"*ANVDNN).LT  -  EPS)  GC  TO  90 
ANV=  (ANVUP+ANVD  **N)*-5 

8  C  CCNTINUE 
SC  CCNTINUE 

N  N  L»=  AN  V 

PFINT  *//.EP.NNU.ITR 


M  i 


■vvrtaacMMii 


o  r> 


in  ■' 


If  (  I TR.EO-HAX I  TP )  PRINT  //#••**  WARNING  •**  THE 
1  ANSRtRMAY  NOT  C  ON  V  EH  GF  * 

Ef=EP«CE 

anv-anvup/2;anvcwn=o 


ICC 

CCNT INUE 

1  C5 

CC  TO 

PUNT 

9 

//*  '  DELAY  TIME  TOO 

LARGE  •  »  T  OEL  A  Y 

11C 

GC  TO 

FF1  NT 

99V 

•  //  »  *  illegal  celay 

FUNCTION* 

SS'1 

SlOP 

END 

SUBROUTINE  DIEFtllS  (T*Y»OY) 


c 

C  Tf-lS  SECTION  PROVIDE  A  SET  OF  FIVE  SECOND  OKOER 

C  E  QUA  HONS 

C  *  HE  R  £ 

C  F  0  -  Y  ( 2  )  G  I  =  Y  (  8  ) 

C  F  l  -  Y  (  A  )  G  2=  Y  (  1  C  ) 

F  2  1  Y  (  6  ) 

CIMENSlCN  YdC)rCY(lO) 

CCHMON  EP »  RK»*B»EPANV 
CCMMON  /SL IST/SC*S1»S2 

C  CNN  ON  /BLAl/  C 1 »C2#C 3,C A*C5*C6*C7 »C3»C9, Cl CrCl 1*C12 
1 »  Cl 5  *  C 1 A»  C15#CI6#C1? 

Z  1  (2 )=  Y(1  ) 

CY(l)  =  -SO*SC*YC2)-EP«(Cl«Y<2)‘Ytl)*C2*(Y(A)*YC3)*Y(S) 

1  «Y(  n  )  ♦  C3*(Y(10)*Y(9)*Y(6)«Y(S)>> 

Cmi--Y(3) 

CYC 3)*-S1*S1*YCA)-EP*<CA«YC2)*Y(3)*C5*YC4)*Y(1)  ♦ 

1  Cb«(Y(8)*Y(9)«Y(A  )•  Y  (5  >>*C7*(Y(b>*V(3>*V(iO>*Y(M>> 
CWb)M(S) 

CY(S)»-S2*S2*Yl6)“EP*(C8*Y(2>*YCS)*C9»Y(b)*Y<l)«  CIO* 
l(  Y<  8  )  *  Y<  2)  -Y(  A  )  •  YC  3))  ) 

[us  >-Y<n 

C  Y  C  F  )  *  “  S 1  *  S  1  *  Y  (  0  )“EP*(C4»Y(2)*Y(7  )  ♦  C  S  *  Y  (  0  )  •  Y  (  1  )  ♦  Cb* 
1(Y(4>«Y(9)-Y(6)*Y(5))*CF*(YCIC>*Y(3)-Y<6)«Y(2>>> 

CYC  13) =Y(9  ) 

C  Y  (  V  )  -  “  S  2*  52*Y(  10)“EP*(Cd*Y(2  )*V(9)*C9*Y<  10  )  *  Y  <  1  )  CIO 
l  »  (  YC  d)* Y<  3  > ♦ Y C A  > • V  C  7))  ) 

IF  (RN.EQ.  2)  RETURN 

C 

C  IFIS  PROVIDE  THE  COMBUSTION  TERMS 

C 

CYC  1  )  s  0  YC l  )aMB*(Y(l  )*EPANV*(CllMlc)*Y(2)*Cl2*(Y(4)*Y 
1(A)*  Y(8)*Y(h))«  Cl 3*( Y (l 0  )* Y  A 1 0  >  * YCb >* Y( b > > > > 

LY(3)  DY(3)-«sj*(Y(3)*EPANV*CClA*Y(2)*YCA)*C15*(Y(A)*Y 
1(f)*  YC  8 ) * Y<  1  C  )  )  ) ) 

CYCS):0Y(‘3)-Wb*CY(,j)«EPANV«(C16*(YlA)*Y(A)-Y(8)*Y(a)) 

1  *C  W*YC  2  )  *Y  (6)  )  ) 


M 


JI^«W  -y  n.  -  • 


o  o  o 
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C1(7)=0Y(7)-W8*(V(7)«EPANV*(C14*Y(2)*Y(6)4C15*<Y(4)*Y 
1(10)  Y(8)*Y(6))>) 

CY(9)*DY<9)-W8*(Y<9)4EPANV*(C17*Y(2)*Y(10)42*C1&*Y(4> 

10(8  ))  ) 

*  E  TORN 
BO 

S13R0LTINE  RKDELY(N,T,Y»DY,H,G»K»EPS»  *»  HHALF  ,K  1 ) 

C 

This  section  perform  a  fourth  orcer  runge-kutta  methgd 

MTH  TIME  DELAY  FUNCTION 

Cl  ME  NS  ION  Y(1Q),CY(10),Y2(10),Y3(10)*G( 10  *  500) 

COMMON  EP»  RK,h6,EPANV 

COMMON  /3LK1/  C1,C2»C3,C4»C5,C6,C7,C8,C9»C1C,C11»C12 
1  »  C  1 3  »  C  1 4  >  C15,C16,C17 

CALL  C I FFU  N ( T , Y , C  Y  ) 

CC  10  1=1, N 

Y2(X)=Y(I  )  ♦  HHALF*(DY(I)  *  G(I,1)> 

G(I,1)=(G( I,l)«G(I,2))/2. 

1C  CCNTINLE 
T  =  T  4Hh  ALF 

CALL  CIFFUN  (T, Y2,DY) 

CC  20  1=1, N  } 

Y  3 ( I  )=  Y  (I  )  ♦  HHALF*(DY(I  )  ♦  G(I,1)) 

2C  Y  2  (  I  )  =  Y2  (  I  )  *  2.  *  Y  3(  I  ) 

CALL  D I FFU  N (  I,Y3,DY) 

CC  40  1=1, N 

Y3(I)=Y(I)  ♦  K*(DY(I)  ♦  G ( 1,1)3 
IF(K.LT.l)  GO  TO  40 
CC  30  J  =  1,K 

3  C  G ( I »  J )  =  G ( I ,j«l) 

4  C  Y2( I )  =  Y2(  l  )  ♦  Y3( I ) 

T  =  T  < HHALF 

CALL  CIFFUN  (T  ,Y3,0Y) 

CC  50  1=1, N 

Y  C I )  =  (Y2(I)  -  Y ( I )  ♦  HhALF*(CY(I)  ♦  G(I,l)))/3. 

5  C  CONTINUE 

ENTRY  DLYFLN  (  N  ,  Y  ,  G  ,  N,  E  P  S  ,K  1  ) 

C 

C  TUS  GIVES  THE  FORM  OF  DELAY  FUNCTION  G 

C 

IF  (K.LT.l  )  RETURN 
IF  (RN.EQ.2)  RETURN 

G( 1 , K 1  )=  W8*  EPANV»(C11*Y(2)*Y(2)4C12*(Y(4  )*Y 

1(A)*  Y(8)*Y(8))«  Cl  3*(  Y(1  C  )*  Y(  10)4  Y(6  )*Y(6)  )  ) 

G  (  3  »K  1  )  =  mE-  EPANV*(C14* Y(2  )* Y(4 )*C15*( Y( 4 )*Y 

1  (£  )♦  Y ( 8  )*Y(1C  ))  ) 

G(5,Kl)=wa*EPANV»(C16*(Y(4)4Y(4)-V(8)*Y(8))4C17*Y(2)*Y 

1(E)) 

G  (7 , K  1  )=  W3*  EPANV*(C14*Y(2)*Y(8)4C15*( Y(4)*Y 

1(10)  Y(8)*Y(6)  )  ) 


■  •  wvrirnmKm*  *4 


G(9*K1)=W3*£PANV*(C17*Y(2)*Y(1O)*2*C16*Y<4)*Y(0)) 

SETURh 

£KD 

SLOCK  DATA 

CCMMON  /0LK1/C 1*CZ*C3  »C4  *C5*C6*C7*C8*  C9*C  10  *C 11 ,C 12 
1»C13*C14*  C1S,C16*C17 

CCMMON  /SL IST/SO  *S  1 *S  2 
CATA  S0*S1, $2/ 3. 03171, 1.04110, 3. Ob  4 24/ 

C  A  T  A  Cl*C2»C3»C4/4. 1373*1. 0423*-. 2034,-1. 9394/ 

CATA  C  5*  C6  »  C  7 *Cfi*C9  / - 2 . 3 1 2  3*  1 - 7 10 / *1 . 4 02 0 , -2 . 705* 

1" 1 .0 108/ 

CATA  C10*C 1I*C12,C13/1.1318,2.506,.48C,-. 196/ 

CATA  C14, C15*C 16 *C17/*2. 4243*1. 3534 44 5*-. 44 7,-3. 481/ 
EM) 


APPENDIX  C 


HOG RAM  FOR  MODAL  AMPLITUDE  AND  PRESSURE 


IQ  6 


....  ,  ■ 


r>  o  o  o  o  (i  o  n  cl  o  o  o  r )  ri  r>  o  o  o  o  r>  n  o  o  o  o  r?  n  r<  r»  o  o  n  o  o  o  <"> 
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iSET  LlSTUTOBlNJ 
i  8 INC  -  F9CH  *F0RTlI8/=  FREE 

FI  IF  1(KInC=0ISK  »MAXRECSlZE  =  15»BLCJCKSlZ£=420»  AREAS  =  4 
1/AREASI2E=1400»FILETYPE=7> 
r 

C  KIN-WING  wONG  JUN301976 

C 

C  CCM8UST10N  INST  ABILITY 

METHOD:  THIS  PRGGRAM  USED  THE  GALERKIN  METHOD  TO 

GENERATE  FIVE  SECOND  ORDER  DIFFERENTIAL  EQUATIONS. 

A  FOURTH  ORDER  RUNGE-KUITA  METHOD  IS  USED  TO  SOLVE 

the  resulting  equations. 

THE  PROGRAM  GIVES  T  HF  MODAL  AMPLITUDE  AND  THE  PER  T  USE  2 
PRESSURE  IN  GRAPHICAL  FORM. 

THERE  AWE  TWO  TYPES  OF  INPUT  QATAS.  THE  FIRST  SET  OF 
INPUT  DATA  USED  FORMAT  FREE  INPUT.  A  COMMA  SERVES 
AS  A  DELIMITER  AND  THE  DATA  ENTER  AS  FOLLOw: 

Tl  INITIAL  TIME 

TF  FINAL  TIME 

DT  STEP  SIZE  FOP  TIME 

NP  PRINT  FREQUENCY 

NDF  1  FGP  THE  DELAY  TIME  APPROCHING  0 
0  OTHERWISE 

IGAS  1  FOR  GAS-DYNAMIC  TERMS 
0  OTHERWISE 

F  0  »  INITIAL  CONDITION  MODIFIER 
OF  EITHER  1  OR  0 

R  RADIUS  OF  THE  CTLINOER 

THE  SECOND  SET  OF  DATA  USED  THE  NAMELIST  (LIST) 

INPUT  OUTPUT  OPTION.  THE  DATA  CAN  8£  ENTEREO  IN 
ANY  ORDER.  HCwlVER  THE  FOLLOWING  DATA  MUST  8E 
PROVIDED  INITIALLY. 

TOELAY  DELAY  TIME  IN  MULTIPLE  OF  PI 

EP  ORDER  PARMETEW 

w8  STEADY  BUPNING  RATE 

ANV  INTERACTION  PARMETER 

PLOT  T  FOR  PRESSURE  PLOT 

F  OTHERWISE 

IF  PLOT  IS  TRUE  /  ONE  MORE  DATA  (NUMPT)  IS  NEEOED  TO 
PROVIDE  THE  NUMBER  OF  POINT  WITHIN  0  AND  2*PI. 

CIm£NSI0N  XX(1CO)/YY(1CO)/DY(10)»P(5)/EO(5)#DF(5)»G 
1C  10/ 5  C  0 ) »  Y (  10) 

CCMMON  EP/IGAS/  RK/W3/EPANV 
CCM-ON  /SL IST/SC/S1/S2 

LC.MMON  /8L  K  1/  C1/C2/C3/C4/C5/C6/C7/C8/C9/C10/C11/C1? 


10*? 


1.C1 3.C14.  C15.C1G.C17 

CCMMON  /BLK 2/C  0  2.  CO  3.  CO 4  .CO  5.  C06.C0  7, Cl 2P.C1 3P.  Cl 4P 
1  *  Cl 5P  »  C22.C23.C24.C25.GAMA 

COMMON  /BLK3/  3J0.BJ1.BJ2.COSC.C0S20.  SIN0.SIN20 
CAT*  PI/3.14159/ 

LOGICAL  PLOT 

NAMELIST  /LIST/  TDEL  A  Y  .  Z  P  .  W  B  *  A  N  9  .  PL  OT 
RE  AD (5./)  T  I  »  T  F  .DT.NP.NOF.IGAS.FO.DF.R 
PRINT  *//.TI.TF.DT.NP.NDF 

IF  (IGAS.EQ.O)  PRINT  //.‘GAS  CYNAHIC  TERM  IS  OFF* 

PRINT  •// »  F 0  »DF  . R 

EFS=1.E-S 

CALL  iiESJ(  SO*R. 0.3J0.EPS. IER) 

CALL  BE  S J  ( S  l  «R  » 1.BJ1.EPS.1ER) 

CALL  BESJ(  S?*R» 2.8J2.EPS. IER) 

999  READ(5»LI ST  »END=99) 

REwINC  l 
TCLY  =  P  I *  TO  E  L  A  Y 
EFAN9 -EP*ANV 

IF  (EPANV.EC.O)  EPANV  = ANv 
*RT TE(b.Ll ST  ) 

N  - 1  0 

hhalf-ot*.5 

K=TDLY/DT 

IF  (M  .GT.  500  )  GO  TO  500 
IF  ( K.LT.l .AND. NDF.EQ.  1  >  K  =  1 
DC  30  I  =  1 »  K 
CC  30  J=  1  »  N 
30  G(J»I>=0. 

DC  10  1*1. N 
1C  Y  C I  ) =C - 
C  INITIAL  CONDITION 

C 

T-TI 

Y(2)=F0(1)»Y(4)=F0(2)1Y(6)=F0(3)AY(3)-F0(4)AY(10)=F0 

1(5) 

Y(1  >  =0F(  1  >  «S0;  Y(  3>=DF(  2)*S  WY  (5>  =  CF  (  3  )*S2 
Y(F>=DF(4>*SI;Y(9)=DF<5>  *S2 
NSTEP=(TF-TI  )/DT 
IF  (.NOT. PLOT)  GO  TO  40 
CALL  PREGUE(Y.P) 

«RI TE(  1  )  T  *  P 
4  C  CONTINUE 

PRINT  //."FUNCTION" 

wRI  TE(6»100)  T.( V < I  ).I  =  2.N.2) 

*.5  CCNTINUE 

CALL  OLYFGN  ( N . Y  .G » K . E PS  * K l  ) 

CC  20  1*1 . NSTEP 

CALL  RK  DEL  Y (N»  T  »  Y.DY.D  T.G.K.EP.  4 40  0.  HH  A L F  .K 1  > 

DC  55  J=2. N*2 


109 


IF(A3S(Y(J  > ).LT •  2  )  GO  TO  55 

A7-2 

55  CONTINUE 

IF  (MOOT  I»NP  ).NE.O  )  GO  TO  20 
IF  (.NOT. PLOT)  GO  TO  60 
CALL  PHESU£( Y*P ) 

WRITE(l)  T .P 

60  nF.I  TE(6#100)  T»(Y(J),J=2#N»2) 

65  CONTINUE 
2  C  CONTINUE 

IF  (RK.GT.O)  PRINT  //»• UNSTABLE* 

IF  (.NOT. PLOT)  GO  TO  70 
ENDF1LE  1 
RERINC  1 

PRINT  //,' PRESSURE  * 

REA0(5»/)  NUHPT 
0  £  =  ° I  *  2./NUMPT 
NLM=NLrtPT*  l 
TE-t°  =  C;XXC  l  )=0I 
CC  7 5  I=2>NUM 
T£MP  =  TENP*  03 
1 5  XX(I)=TEMP 
96  REAO(l#ENO  =  rO)  T  » P 
PRINT  //»"TIM£=  " »  T 
SLMY=C. 

CC  80  1=1#  NUN 
3C=XX( I ) 

Y  Y  (  I  )  =  -r.AMA*(P(l)»P(2)*CGS(3C)*P(T)*C0S(2*QQ)*P(4) 
l  *  SI  N ( 20  )♦  P(5)*SIN(2*3C) ) 

SU* Y=SUNY*  Y  Y  (  I  ) 
eO  CONTINUE 

IF  (SUMY.NE.O.)  CALL  PL0T2D(XX>YY»NUM*1Q0/100) 

GC  TO  98 
70  CONTINUE 
R  I*  =  0 

GC  TO  999 

<•00  PRINT  *//»  *  ILLEGAL  OELAY  FUNCTION* 

GC  TO  99 

5  CC  PRINT  * // # • OEL  A  Y  TIME  TOO  LARGE* 

99  S10P 

100  F  CRN  A  T  (X.8E12.5) 

11C  FORMAT ( 13X  .7E12.5) 

E  NO 

SUBROUTINE  DIFFUN  (T,Y*OY) 

C1MENSI0N  Y(1G>#0Y<10) 

CCMMON  EP# I G A  S  »  RN»W9»EPANV 
CCMMON  /SL IST/SG»S1*S2 

CCMMON  /BL  N 1 /  Cl»C2*C3rCA»C5*C6»C7*C8»C9»C10*Cll#C12 
l»Cl  3  »  C  1  <• »  C15»C16*C17 

C 

CY(2 )  =  r < 1  ) 


C  Y<  1  )=  “SO*SO*  Y  (2  >-EP*IGAS*(Cl»Y<2)*Y(  1>*C2*(Y(<.)*Y<  3) 
1«Y(3>*Y(7>>  •C3*(Y(10)«Y(9)«Y(6)*Y(5>>) 

DY(<*)  =  Y(  I) 

CY(3  )  =  “Sl*al*Y(4  )_EP*IGAS*(C4»Y(2)*Y(  3)*Cc'*Y(4)*Y(  1  )  ♦ 

1  C6«<Y(8>*Y<9)*YU)*Y(5>>*C7*CY(6)*V(3>*Y110)*Y 

2  C  7  >  )  ) 

2  Y<S  >=  Y<5  ) 

0  Y  (  5 ) =”S2*$2*Y£  6 )“E  P*  I GAS»(C8*Y( 2  )  *  Y(  5>*C9*Y(  6)*Y(  1)* 

1  C10*( Y(8)*YC7)-YC4)*Y( 3  ))) 

0Y(8)=Y(7) 

0  Y  (  7)  =  ~S1*S  1  *Y  <  8  )"FP*I  GA  S  •(  C4  *Y(  2  )•  Y (  7  >*C5*Y(3  )*Y(1  )♦ 

1  Cb*(Y(4)»Y(9)-Y(8)*Y(5))4C7*(Y<10)*Y(3)-Y(6)*Y<7))) 
CY( 10  )  =  Y  (  9  ) 

0Y( 9 )=“S2*:>2*Y(  10)~EP*IGAS*(C6*Y(2 ) * Y C9)*C9*Y( 10  )  *  Y  (  1  ) 

I  Cl  0  •  (  Y(  6  )  *  Y(  3  )♦  Y<  «  )•  Y(  7  )  )  ) 

IF  (RIS.EQ.  2)  RETURN 

CY(l>-0Y(l)-rfS*(Y(l)*CPANV»(Cll*Y(2)*Y(?>*C.  12*( Y( *  )*Y 

I I  ♦  Y(8)*Y(8)>*  Cl 3*< Y<1 0  )* Y(1 0 ) »Y(6  )* Y(6 > ) > ) 

CY<3)  =  DY(3)-Rd*<Y(3MEPANV*<cn*Y(2)*Y(4)*C15«(Y(<.>*Y 

1  C  £  )♦  Y( 8  )  * Y(  l  0 ) ) )  1 

CY<5)=0Y(5)-W3*(Y(5)«E°ANY*<Clb*(Y(4)*Y(4)-Y(3)*Y(8)> 
l  *C1 7* Y(2  )*Y(b)  )) 

3Y(7)=0Y(7>-wd*(Y(/Y*roANV*(Cl‘*»Y(2>*YC5>-*Clb*(Y(A)*Y 
1<  10)  Y( 8) *Y(fc ) >)> 

CY(9  )-OY(9  )“WB*  (  Y(9  )»Ep>ANV*tC  l  7  *  Y  (  2  )  •  Y  (  10)*  2  •  C  1  b  *  Y  (  4 ) 

1  *  Y  t  8  >  >  > 

RETURN 

ENJ 

SUBROUTINE  RKl)ELY<N,T,  Y,DY  ,H,G*iU£PS*  *  ,  H  Ha  L  F  >K  1  ) 

C  HE  NS  ION  Y(10)»OY(10)»V2(10)»Y3£10)»G£  10  »bOO ) 

CCMMQN  EP  » 1  G  A  S  ,  RK»wB»EPANV 

CCHMON  /8LKI/  Cl ,C2»C 3»C4,C5, C6,C/,C8,C9, C10,C 1 1 »C 12 
1*C13,CU»  Clr>»Clb,C17 
CALL  OIFFGN<T*Y,DY) 

3C  1  1=1, N 

Y2(1)  =  Y(I)  ♦  HHALE*(0Y(1)  *  G( 1 ,  1 )  ) 

3(1 , 1 ) = ( G( I,l)«G(I,2))/2 
CCNT INUE 
T  =  T  ♦ HH ALF 

CALL  0IFFUN(T»Y2,DY> 

00  2  1=1, N 

Y3U)=Y(I>  ♦  HHALFMDY(l)  ♦  G(  1  ,  1  >  > 

Y21  I  }=  Y C  I  )  ♦  2*Y  3<  I  ) 
call  CIFFUNt  T,YJ»DY) 

OC  1  1=1, N 

Y  3(  I  )  =  Y  (  I  )  ♦  H  *  (  0  Y  (  I  )  ♦  GC  1 , 1  >  > 

IF(N.LT.l)  GO  TO  3 
OC  10  J= 1 , K 
G  (  I  »  J  )  =  &{  I  ,  J»1  ) 

Y2(I  )  =  Y 2 ( I )  ♦  Y  3 ( 1 ) 

T  =  T«  HHALF 


CALL  CiFFUN  (T  .Y3.0Y) 

CC  4  I  -  1  *  N 

Y(I)  =  l  Y2(  I )  -  Y { 1 )  ♦  HHALF*(DY(I)  ♦  G(I,l))>/3. 
CCNT INUE 

►NTRY  GLYFUN  (  N.Y.G.K.EPS.K1  ) 

IF  (  K  .  L  T  .  1  )  RETURN 
IF  ( RK.EQ. 2)  RETURN 

G(1»KI>  =  *8*  EPANV*(C  ll*Y(2  )*Y(2  >  ♦  C  1  2  *  C  Y(  4  )*Y 

1(A)*  Y<«?)  *Y(8)  >♦  C13*(Y(lC)*Y(lO)«Y(b)*Y(b>)) 

3(3. Nl)*  W8*  EPANV*(C  1  4  *  Y  (  2  )•  Y  (  4  )*C15*(Y(4  )  •  Y 

1(E)*  Y( 3 ) * Y( 1 0  )  )  ) 

G(5.Kl)-Mb.EPANV*(C16*(Y(4)*Y(4)-Y(8)*Y(8))*C17*Y(2)»Y 
l  C  fc  )  ) 

3(7»K1)=  WO*  EPANV*(C14M(2)*Y(8)«C15*(Y(4)*Y 

1(10)  Y ( 8 )  • Y ( b  )  )  ) 

G(9,M)=W3*EPANV*(Cl7*Y(2)*Y(10)*2*Clb*Y<4)*Y(8>) 

RETURN 

END 

SUBROUTINE  PRESUE(Y.P) 

THIS  SECTION  PEKFGRM  THE  PRESSURE  CALCULATION 

DIMENSION  Y(10)»PC5) 

COMMON  EP  »  I  GAS.  RK»*n.EpANY 

COMMON  /8LN2/CQ2»C03»C04.C05»COb»C07*C12P»C13P»C14P 
l.C15p.  C22.C2  3.C2-*  .C25.GAMA 

COMMON  /BL  N  3/8  JO.  8  Jl  .  B  J2.C0SJ  .C0S20  .SI  NO.  SI  N2Q 
0  F  0  -  Y  (  1)»L)F1=Y(  3)pr)F2=Y(S);DGl=Y(7);0G2=Y(9) 

F0^Y(2);  F  I  -  Y  (  4  )/  F?-Y(6)i  CI  =  Y(c)J  G2=Y<10) 
P(l)=riJ0*(DFO*EP*(CO2*FO«FO*C03*(Fl*Fl*Gl*Gl)«CO4*(F2 
1*F2*G2*G2)  ♦  C05*DFO*DF0*CO6*(CFl*i)Fl*t)Gl*0Gl)*C07* 

7(lF?*CF2*  DG  2  *DG?  )  )  ) 

P(2)=3JI*(0F1*EP*(C12P*FC*F1*C13P*(F1*F2*G1*G2)*C14P 
!*0Fa*3Fl*  C15P*(0rl»0F2*0Gl*DG2))  ) 

P(3)=BJ2*(0F2»EP*(C22*F0*F2*C23*(FI*F1-G1*G1)«C2A*0F0 
1 *CF  2  ♦  C25«( GF1*DF 1-0G1*0G1  > >  ) 

P(4)  =  eJl*(DGl*£P*(CI2P»F0*Gl«C13P*(Fl*G2“Gl*F2)*Cl4P 
1*CF0*0G1»  C15P*(0F1*0G2-CF2»CG1))) 

p  (5  ) -6 J2*<  OG2*EP«(C?2*FO«G2*2 . *C23*F1 *G1*  C24*OFO*OG2* 

1  2  •  *  C2  5  *  OF 1 »OGI  )  ) 

RETURN 

ERO 

3  L  J  C  K  OATA 

COMMON  /BL  K  1  /  C 1 .C2.C 3. C4.C5.CS.C7.C8. C9.C10. Cl 1.C12 
1.C1 5.C14,  C15.C1S.C17 

COMMON  /SL 1  ST/SO  , SI »S2 

COMMON  /'3LK2/C02.CG3»C04»C05»C0b»C07»C12P»ClSP.Cl4P 
1.C15P.  C22.C23.C24  .C25.GAM  A 

CAT  A  S0.S1 »  S2/ 3.83171.1.84113.3.05424/ 

GAT A  Cl.  C2.C3.C4/4. 1373. 1.0423.-. 2084. -1.9394/ 

CATA  C5.C6.C7.C3.C9  / - 2.  31  2 3 .  1 . 7 1 8 7  . 1  . 4  82 8 . “2 . 7 35 . 


1*3.0568/ 

Data  C  1 0  *  C  1  l  »  C  1  ?  ,C13/1.1  318,7.586,. 463  >»“•  196/ 

:*TA  C14,C15,C16,C1  7/-?. 4243, 1.8534445*-. 447,-3. 481/ 
CATA  C0  2, C  0  3»CC4»C05»rC6»C0//  I  . 29  50  *0.2400, -.0982, 
l-.irb2  ». Ou 07, . 0494/Cl 2P*C13P»C14P,C15P/ -1.2121 

2»  .9267 ,.1651,“. 1 98 7/ C 22 , C2 3 »C 2 4 , C25 /- 1 . 7 4  Ob , -.  22 35 
3,  .237  1  ,-.l  754./ 

CATA  CAMA/1 .2/ 

ENC 


